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ABSTRACT 


In  this  project,  computational  tools  were  developed  that  help  in  designing  and 
analyzing  multi-functional  composite  structures  that  have  sensing  and  actuation 
capabilities.  Magnetic  actuation  of  composite  structures  using  built-in  or  embedded 
electromagnetic  devices  was  studied  and  a  method  for  detecting  load  and  damage  in 
composite  structures  by  determining  change  in  resistivity/conductivity  was  studied.  An 
algorithm  for  solving  inverse  problem  to  determine  average  resistivity  values  in 
composite  structures  was  demonstrated.  The  method  can  use  data  from  an  arbitrarily 
large  number  of  electrodes  to  compute  average  values  of  resistivity  or  conductivity  for 
the  structure.  Finite  element  models  for  the  structure  are  used  to  solve  the  forward 
problem,  making  this  method  very  general  and  applicable  to  arbitrary  shaped 
structures.  Ideally,  the  electrodes  should  be  embedded  in  the  structures  during  the 
manufacturing  process  itself  so  that  it  can  be  used  for  quality  control,  detection  of 
defects  as  well  as  subsequent  health  monitoring. 

The  idea  of  using  magnetic  forces  to  actuate  structural  mechanisms  was  studied. 
The  main  application  of  interest  is  micro  air  vehicle  wings  that  are  shell  like  structures. 
Topology  optimization  method  was  studied  as  a  potential  method  for  designing 
structures  that  have  specified  modes  of  deformation.  The  structure  is  then  actuated 
using  magnetic  actuation  means  built  into  or  around  the  structures.  Several  actuators 
were  studied  including  solenoid  actuators  and  coil  actuators.  After  systematic 
comparison  of  several  designs,  it  was  concluded  that  a  coil  actuator  built  into  composite 
structures  is  an  ideal  means  for  actuation  of  composite  structures.  A  conceptual  design 
of  a  flapping  wing  vehicle  was  developed  that  is  designed  to  actuate  by  built  in 
actuation  capability  of  the  body,  wing  and  support  structures.  No  external  mechanisms, 
motors  or  linkages  are  needed. 

Computational  tools  were  developed  to  design  and  analyze  structures  actuated  by 
magnetic  forces.  Magnetostatic  analysis  capability  was  implemented  into  a  pre-existing 
software  (named  IBFEM)  developed  at  the  University  of  Florida  that  can  perform  finite 
element  analysis  without  the  need  for  generating  mesh.  Solid  and  surface  geometry 
modeled  on  commercial  CAD  software  can  be  imported  into  this  software  and  analysis 
can  be  performed  without  approximating  the  geometry  using  a  conforming  mesh.  The 
structured  mesh  approach  has  been  demonstrated  to  work  for  magnetostatic  analysis 
and  validated  using  several  examples  with  known  solutions.  The  approach  has  been 
demonstrated  for  both  2D  and  3D  magnetostatic  models.  Structured  mesh  is  easy  to 
generate  and  the  elements  are  regular  and  not  distorted  as  in  traditional  finite  element 
mesh.  Magnetic  forces  were  computed  by  integrating  the  magnetic  force  density.  These 


1 


forces  are  then  used  in  a  subsequent  structural  analysis  to  determine  the  deflection  of 
the  structure.  Shell  elements  based  on  uniform  B-spline  shape  function  were 
implemented  into  IBFEM.  One  of  the  key  advantages  of  using  these  elements  is  that  a 
structured  mesh,  which  is  easy  generate  automatically,  can  be  used  for  the  analysis. 
Both  quadratic  and  cubic  B-spline  shape  functions  were  tested  and  it  was  found  that 
cubic  elements  provided  very  good  results  with  fewer  elements  than  quadratic 
elements.  Computational  cost  is  higher  for  these  elements  compared  to  traditional  shell 
elements  but  often  fewer  elements  are  needed  to  get  accurate  results  with  cubic 
elements.  The  time  taken  to  create  the  model  is  significantly  lower  because  structured 
mesh  generation  is  easily  automated. 
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I.  INTRODUCTION 


Composite  materials  have  become  the  structural  material  of  choice  in  many 
aerospace,  automotive  and  other  applications  where  low  weight,  high  strength  and 
rigidity  are  required.  For  applications  such  as  in  micro  and  unmanned  air  vehicles,  in 
addition  to  the  above  requirement,  it  is  necessary  to  pack  more  and  more  functionality 
into  less  and  less  space.  This  provides  the  motivation  to  make  the  structure  of  the  air 
vehicles  multifunctional,  allowing  it  to  perform  important  tasks  including  actuation, 
sensing,  energy  storage  and  energy  harvesting  in  addition  to  providing  structural 
support  and  rigidity. 

Multifunctional  composites  can  significantly  increase  the  duration  of  flight  and  the 
payload  of  unmanned  air  vehicles  (UAVs),  especially  micro-air  vehicles  (MAVs).  Most 
of  the  research  related  to  multifunctional  materials,  in  the  past,  has  been  in  the  area  of 
specialized  materials  such  as  piezoelectric  and  magnetostrictive  materials  that  have 
been  used  to  design  and  build  both  actuators  and  sensors.  These  specialized  materials 
are  typically  not  good  structural  materials  and  therefore  are  only  used  as  attachments 
on  a  structure  to  achieve  sensing  capability  or  to  induce/damp  vibrations  etc.  Effect  of 
large  currents  through  metallic  plates  and  shells  have  been  studied  in  the  past  to 
explore  potential  applications  such  as  providing  better  impact  resistance  to  armor 
plates.  Conductors  carrying  large  currents  through  a  magnetic  field  experience  a 
damping  force  that  could  potentially  provide  improved  impact  resistance.  However, 
metals  lose  strength  at  higher  temperatures.  Therefore  the  amount  of  current  that  can  be 
applied  is  limited  by  Joule  heating  of  the  plates.  For  composite  materials  also  heating  is 
a  problem  that  limits  the  current  that  it  can  carry  so  that  electromagnetic  damping  may 
not  provide  significant  impact  resistance.  However,  heat  generated  by  currents  can  be 
used  beneficially  to  provide  better  curing  of  composites  as  well  as  for  self  healing  after 
an  impact  especially  for  composites  that  have  thermoplastic  polymer  matrix. 
Experimental  evidence  suggests  that  electric  currents  passing  through  composite 
structures  improve  their  impact  resistance  [1-2]. 

In  this  project,  methods  for  designing  and  fabricating  multifunctional  structures 
were  studied.  In  particular,  the  focus  was  on  designing  sensors  and  actuators  that  can 
be  of  great  value  to  applications  such  as  micro  air  vehicles.  Feasibility  of  sensing 
applied  loads  and  damage  on  composites  was  studied.  The  idea  of  using  embedded 
circuits  and  magnetic  forces  for  actuation  of  structures  was  also  studied.  Analysis  tools 
and  design  methods  were  developed  for  designing  structural  actuators  that  are 
structural  mechanisms  actuated  using  magnetic  forces  generated  by  embedded  circuits 
and  magnets.  The  main  activities  are  summarized  below. 
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1.  Sensors 


The  conductivity  /  resistivity  of  composite  plates  were  studied  with  the  goal  of 
measuring  changes  in  these  properties  as  a  means  of  detecting  applied  loads  and  for 
sensing  damages  within  the  structures.  For  composite  plates  and  shells,  the  most 
common  types  of  damage  are  due  to  delamination  which  occurs  between  layers  and 
fiber  breakages.  Both  these  types  of  damage  may  not  be  visually  detectable  from 
external  appearance.  An  inverse  problem  solving  algorithm  was  developed  for 
determining  the  overall  resistivity  of  a  composite  structure  using  voltages  measured  at 
external  electrodes.  This  provides  a  means  of  detecting  damage  by  comparison  with 
resistivity  values  of  undamaged  material.  The  measurement  accuracy  was  sufficient  to 
detect  damage  due  to  delamination  and  fiber  breakage.  Strains  also  cause  small  changes 
in  resistivity.  The  inverse  method  was  found  to  be  not  accurate  enough  to  determine  the 
magnitude  or  location  of  the  applied  load.  The  main  problem  is  that  the  resistivity 
changes  due  to  applied  loads  are  very  small  and  yet  highly  nonlinear  and  lacking 
repeatability. 


2.  Actuators 

Magnetic  fields  can  be  created  by  currents  flowing  through  fiber  networks  or 
embedded  conductors  within  composites.  These  magnetic  fields  can  be  used  to  actuate 
or  deform  the  structure  which  in  turn  functions  as  a  mechanism.  This  method  can 
provide  means  to  create  very  compact  devices  that  have  the  capability  to  deform  or 
serve  as  actuators.  The  potential  applications  include  morphing  (or  shape  change)  of  air 
vehicles  and  air  vehicles  that  fly  by  flapping  their  wings.  A  key  challenge  in  designing 
such  structures  is  the  lack  of  suitable  design  and  analysis  tools.  Many  commercial 
programs  provide  electrostatic,  magnetostatic  and  electrodynamic  analysis  capability. 
However,  the  ability  to  perform  coupled  magneto-elastostatic  analysis  is  needed  for  this 
application.  Furthermore,  the  geometry  of  flapping  wing  type  structures  are  complex 
requiring  shell  like  analysis  for  the  composite  wings  while  using  solid  structures  to 
model  the  magnetic  actuators.  Creating  models  for  such  systems  is  difficult  in 
conventional  finite  element  analysis  software.  A  method  that  does  not  require 
conforming  mesh  was  developed  for  coupled  magneto-elastostatic  analysis  of  multi¬ 
material  systems  by  extending  Implicit  Boundary  Finite  Element  Method  (IBFEM).  This 
method  of  analysis  was  developed  at  the  University  of  Florida  and  has  been  used  for 
linear  elastic  and  heat  transfer  analysis  is  the  past.  The  method  was  extended  to  enable 
magneto-elastostatic  analysis  as  part  of  this  research  and  applied  to  study  several 
actuator  designs  for  flapping  wing  design. 
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3.  Coupled  Magneto-elastostatics  using  IBFEM 

Implicit  boundary  finite  element  method  (IBFEM)  avoids  the  need  for  generating 
conforming  mesh  by  using  structured  mesh  for  the  analysis.  A  structured  mesh  (also 
referred  to  here  as  a  grid)  is  a  non-conforming  mesh  that  is  made  up  of  regular  shaped 
elements  (rectangles  for  2D  or  cuboids  for  3D)  and  is  easier  to  generate  than  traditional 
finite  element  mesh.  The  geometry  of  the  analysis  domain  is  represented  using 
equations  that  are  independent  of  the  grid.  Boundary  conditions  are  applied  using 
solution  structures  that  are  constructed  using  approximate  step  functions  of  the 
boundary  such  that  these  boundary  conditions  are  guaranteed  to  be  enforced.  A  variety 
of  interpolation  functions  and  approximations  such  as  B-splines  can  be  used  with  this 
approach.  IBFEM  was  extended  to  perform  magnetostatic  analysis  and  to  compute 
forces  due  to  magnetic  field.  Using  these  forces  in  a  subsequent  elastostatic  analysis,  it  is 
possible  to  simulate  the  deformation  produced  by  the  structure. 

4.  Analysis  of  shell-like  structures  using  iBFEM 

Shell-like  structures  are  modeled  in  traditional  finite  element  method  using  shell 
elements.  The  geometry  for  such  structures  is  modeled  using  surfaces  that  represent  the 
mid-plane.  In  order  to  avoid  mesh  generation  on  a  surface,  the  Implicit  Boundary  Finite 
Element  Method  (IBFEM)  was  extended  for  the  analysis  of  shell-like  structures.  Three 
dimensional  elements  that  use  uniform  B-spline  approximation  schemes  are  used  to 
represent  the  displacement  field.  The  surfaces  representing  the  shell  passes  through 
these  elements  and  the  equations  of  these  surfaces  are  used  to  represent  the  geometry 
exactly.  B-spline  approximations  can  provide  higher  order  solutions  that  have  tangent 
and  curvature  continuity.  Numerical  examples  are  presented  to  demonstrate  the 
performance  of  shell  elements  using  IBFEM  and  B-spline  approximation.  Models  of 
flapping  wings  were  created  and  analyzed  to  determine  deflections  due  to  magnetic 
forces  produced  by  the  actuators  that  were  designed.  Mode  shapes  of  vibration  or 
oscillation  of  several  wings  designs  were  also  studied  using  this  analysis  tool. 

5.  Design  of  structural  actuators 

Several  magnetic  actuator  designs  were  studied  as  potential  actuation  mechanism 
within  structural  actuators.  Three  traditional  designs  were  studied  first:  Solenoid 
actuator,  clapper  actuator  and  coil  actuators.  Potential  for  using  these  basic  designs  to 
actuate  a  flapping  wing  mechanism  was  studied.  The  coil  actuator  was  found  to  be  the 
most  promising  approach.  Structural  designs  suitable  for  the  flapping  wing  mechanism 
were  explored.  Topology  optimization  technique  appears  to  be  a  promising  tool  to 
design  such  actuation  mechanism.  The  basic  methodology  for  2D  structural  mechanism 
design  using  topology  optimization  using  IBFEM  was  demonstrated.  Further  research  is 
needed  to  extend  these  ideas  to  design  shell-like  composite  structures. 
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II.  COMPOSITE  STRUCTURAL  SENSORS 


1.  Overview 

A  method  to  estimate  the  resistivity  of  composite  structures  using  an  inverse 
problem  solving  algorithm  was  developed  that  uses  voltage  distribution  on  the 
structure  as  data.  Electrodes  attached  to  the  surface  of  the  structure  are  used  to  obtain 
voltage  data  in  response  to  current  injection  through  a  pair  of  these  electrodes.  The 
forward  problem  involves  using  the  finite  element  method  to  predict  the  voltages  at  the 
electrodes  using  known  values  of  resistivity.  The  inverse  problem  involves  solving  for 
the  resistivity  values  using  the  experimentally  measured  voltage  data.  If  the  material 
does  not  have  uniform  properties,  the  computed  resistivity  values  are  average  values. 
Damage  or  defect  in  a  composite  structure  can  significantly  alter  the  average  resistivity 
of  the  structure.  To  explore  the  possibility  of  using  this  approach  to  detect  defects  in 
manufacturing  or  damage  due  to  loading,  the  effect  of  artificially  induced 
damage/defect  on  the  overall  resistivity  of  the  structure  was  studied. 

2.  Sensing  capabilities  of  carbon  fiber  composite  structure 

Carbon  fiber  based  composites  are  of  interest  in  multifunctional  and  smart 
structure  design  because  they  are  conductive  and  there  is  a  correlation  between  changes 
in  electrical  properties  and  applied  strains.  The  resistance  changes  under  a  variety  of 
load  types  including  tension/compression  [3]-[8],  bending  [9],  and  impact  [10]-[16]  have 
been  studied  in  the  past.  The  motivation  for  these  studies  has  been  to  explore  the 
possibility  of  using  carbon  fiber  composite  as  strain  or  stress  sensor  by  measuring  its 
change  in  resistance  due  to  applied  strain.  The  electrical  properties  are  also  affected  by 
any  damage  such  as  delamination  and  cracks  [18]-[27]  that  may  occur  in  the  structure. 
This  provides  a  mechanism  to  sense  damage  by  measuring  the  change  in  resistance  and 
to  identify  delamination  or  crack  or  even  to  quantify  the  energy  of  the  impact  that 
caused  the  damage.  Both  AC  and  DC  measurements  have  been  used  as  means  of  non¬ 
destructive  testing,  damage  detection  and  monitoring  [18]-[19].  Due  to  its  excellent 
mechanical  properties,  carbon  fiber  composites  are  widely  used  as  a  structural  material. 
If  they  can  also  serve  as  a  sensor  that  can  detect  applied  loads  or  internal  damage 
without  the  need  for  external  sensors,  then  they  can  serve  as  a  multifunctional  or  smart 
structural  material. 

The  resistivity  of  carbon  fiber  composite  material  is  orthotropic  and  therefore  it  is 
characterized  by  the  three  principal  values.  In  order  to  measure  the  resistivity  of  carbon 
fiber  composite,  several  specimen  shapes  and  electrode  placements  schemes  have  been 
studied  [3],  [5],  [17].  The  simplest  scheme  would  be  to  apply  a  uniform  current  density 
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on  a  composite  specimen  between  electrodes  on  parallel  faces  and  measure  the  voltage 
at  these  electrodes.  This  is  the  two-probe  approach  where,  current  injection  and  voltage 
measurements  are  made  at  the  same  pair  of  electrodes.  Using  resistance  measured 
between  two  electrodes,  the  resistivity  can  be  calculated  when  the  specimen  has  simple 
geometry  of  known  dimensions  [17]-[21].  Fig.  1  illustrates  both  the  two-probe  method 
and  the  four-probe  method.  The  two-probe  method  is  highly  sensitive  to  the  contact 
resistance  at  the  current  injection  electrodes  because  the  measured  voltage  difference 
includes  the  voltage  drop  across  the  electrode  and  its  interfaces.  In  the  four-probe 
method,  one  pair  of  probes  is  used  for  the  current  injection  at  a  pair  of  electrodes  while 
the  other  pair  is  used  for  voltage  measurement  on  a  different  set  of  electrodes.  Again, 
the  resistance  is  determined  using  Ohm's  law  and  the  resistivity  can  be  calculated  if  the 
applied  current  is  uniform  and  the  specimen  dimensions  are  known.  The  four-probe 
method  provides  results  that  are  more  reliable  because  it  is  not  sensitive  to  contact 
resistance.  These  approaches  for  determining  the  resistivity  are  suited  for  simple  block 
like  specimen  subjected  to  uniform  current  where  the  electrodes  cover  an  entire  side  of 
the  specimen.  It  is  beneficial  to  be  able  to  measure  resistivity  for  arbitrary  shaped 
structures  that  may  be  subjected  to  loads.  Applied  loads  on  the  specimen  can  cause 
contact  degradation  between  electrodes  and  the  carbon  fiber  composite  specimen  [10]- 
[14].  It  is  preferred  that  the  electrodes  can  be  placed  anywhere  on  the  specimen  so  that 
they  can  be  placed  at  locations  least  affected  by  the  loads. 


Fig.  1  Two  conventional  methods  to  measure  a  resistivity:  (a)  Two-probe  method  (b) 

Four-probe  method. 

An  approach  for  estimating  the  average  resistivity  of  arbitrarily  shaped  carbon 
fiber  composite  structures  using  inverse  method  was  developed  as  part  of  this  research. 
An  inverse  method  involves  computing  the  material  properties,  in  this  case  resistivity, 
by  searching  for  values  of  the  properties  such  that  a  model  using  these  values  can 
correctly  predict  a  set  of  experimentally  measured  responses.  Inverse  methods  have 
been  used  in  the  past  to  detect  damage  in  composites.  Todoroki  et  al  [20]-[23]  used  an 
array  of  electrodes  placed  along  the  top  surface  of  a  plate-like  specimen  to  detect 
damaged  areas  using  inverse  method  and  response  surface  models.  The  electric 
potential  computed  tomography  (CT)  approach  [24]-[27]  also  uses  inverse  method  for 
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defect  identification  (delamination  or  surface  crack)  using  passively  observed  electric 
potential  values  on  a  thin  piezoelectric  film  attached  on  the  surface  of  the  structure. 

In  this  project,  an  inverse  problem  solving  method  was  used  to  compute  the 
overall  or  average  resistivity  values  of  arbitrarily  shaped  composite  structures  using 
voltages  measured  on  surface  electrodes  as  data.  If  the  resistivity  of  the  undamaged  and 
defect  free  material  is  known  then  the  average  resistivity  computed  using  this  method 
for  a  given  structure  can  be  compared  with  the  known  values  to  determine  if  the 
structure  is  free  of  defects  or  damages.  The  forward  problem  is  the  finite  element  model 
of  the  composite  structure  that  can  predict  voltages  on  the  set  of  electrodes  distributed 
over  the  specimen  if  the  resistivity  is  known.  The  inverse  problem  involves  solving  for 
the  resistivity  values  using  experimentally  measured  voltages  at  the  electrode.  The 
primary  advantage  of  using  a  finite  element  model,  as  opposed  to  an  analytical  model, 
in  the  forward  problem  is  that  the  specimen  can  then  be  of  arbitrary  shape  giving  us  the 
flexibility  to  determine  average  material  properties  of  real  structural  components  that 
are  in  use  in  automotive  or  aerospace  structures. 

3.  Experimental  Procedures  and  Analysis 

In  this  study,  both  unidirectional  carbon  fiber  composite  plates  and  woven  fiber 
composites  were  used  with  electrodes  attached  on  the  top  and  bottom  surfaces.  Voltage 
differences  between  pairs  of  electrodes  were  measured  using  one  set  of  probes  while 
current  was  injected  at  a  different  pair  of  electrodes  using  a  different  set  of  probes.  Fig. 
2  shows  the  schematic  diagram  of  the  experimental  setup  with  a  specimen  that  has 
eight  electrodes,  with  four  electrodes  on  the  top,  and  the  other  four  electrodes  at  the 
bottom.  As  shown  in  the  figure,  a  current  source  and  voltmeter  are  attached  to  a 
multiplexer  that  can  be  used  to  measure  the  voltage  between  any  two  pair  of  electrodes 
while  current  is  injected  between  any  other  pair  of  electrodes.  Therefore,  for  any  pair  of 
electrodes  between  which  current  is  applied  we  can  measure  voltage  difference 
between  all  the  other  combinations  of  electrodes  available.  Then  by  changing  the 
current  injection  electrode  pair,  even  more  data  can  be  obtained  that  can  be  used  in  the 
inverse  method  to  obtain  an  average  value  of  resistivity  in  the  three  principal  directions. 
This  method  will  be  robust  even  when  the  structure  is  under  any  loading  because  it  is 
not  necessary  to  generate  uniform  current  distribution  between  any  electrode  pairs  and 
electrodes  could  be  located  where  it  is  relatively  safe  from  damage  due  to  any  external 
load. 
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Fig.  2  A  schematic  diagram  of  the  experimental  setup 


As  mentioned  earlier,  the  model  used  for  the  inverse  problem  is  a  finite  element 
model.  A  finite  element  mesh  that  represents  the  geometry  of  the  specimen  with 
reasonable  accuracy  is  needed.  Using  this  FEA  model,  we  solve  the  forward  problem, 
which  involves  finding  the  voltage  distribution  using  the  known  value  of  applied 
current  and  an  estimate  (or  guess)  of  the  resisitivity  values  in  the  three  principal 
directions.  Accurate  values  of  the  resistivity  are  then  calculated  by  minimizing  the  error 
between  the  computed  voltages  at  the  electrodes  and  the  measured  voltages  by  varying 
the  resistivity  values  iteratively.  The  optimization  process  involves  solving  the  forward 
problem  repeatedly  until  good  estimates  of  the  resistivity  values  are  obtained.  Gauss- 
Newton  algorithm  was  used  for  minimizing  the  least  square  error  [28]. 

E»«cUods  2  (ToA)and  6'<9ottom> 

Etactro<i9  1  (Top;  an4  S  (Bottom) 


€l«cirdd«  4  (topland  8  jBotipinl 


il«ctrod«.3(Top)an<>  7(Bottoni) 


Fig.  3  FEA  model  of  specimen  with  eight  electrodes 

Fig.  3  shows  the  FEA  model  used  for  the  eight  electrodes  specimen.  The  size  of 
specimen  modeled  here  is  54x52x1.5  mm^.  The  model  has  3024  hexahedral  elements  so 
that  the  average  element  size  is  2x2x1. 5  mm^.  Electrodes  1  and  7  were  used  as  the 
current  injection  electrode  pair,  with  electrode  7  as  source  and  electrode  1  as  drain  or 
reference  point.  The  current  injection  electrodes  were  modeled  as  surfaces  (or  faces  of 
elements)  that  have  a  specified  normal  component  of  current  density.  The  voltage  at  the 
reference  electrode  (or  associated  nodes)  was  set  to  zero.  For  this  analysis,  the 
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amplitude  of  injected  current  applied  was  28.9mA  and  the  resistivity  was  0.01  m0.m  in 
all  the  three  principle  directions.  Fig.  4  shows  the  voltages  computed  by  solving  the 
forward  problem  using  FEA. 

1 
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Fig.  4  The  voltage  distribution  as  the  solution  of  a  forward  problem 

4.  Inverse  problem 

The  inverse  problem  is  an  optimization  problem  where  the  objective  is  to 
minimize  the  error  between  the  outcome  predicted  by  the  forward  problem  (or  the 
model)  and  the  experimentally  measured  outcome.  For  the  problem  of  interest  in  this 
project,  the  forward  problem  involves  solving  the  governing  equations  of  electrostatics 
using  the  finite  element  method  to  obtain  voltage  distribution  in  a  composite  structure. 
The  governing  equation  and  the  natural  boundary  conditions  are 


V-([o]-V<zi)  =  0  in  Q 

(2.1) 

=  -  J  •  n  =  [o]  •  •  n  on  0Q 

(2.2) 

where,  (j)  is  the  electric  potential  or  voltage,  is  the  normal  component  of  the 
current  density  along  the  boundary  dCl  and  n  is  the  outward  unit  normal  vector  at  the 
boundary,  [o]  is  the  conductivity  matrix  of  the  material  which  is  a  diagonal  matrix 

whose  components  are  the  inverse  of  resistivity  values  in  the  principal  directions.  If  the 
material  is  isotropic,  conductivity  can  be  treated  as  a  scalar  a . 

The  optimization  problem  for  the  inverse  method  is  the  minimization  of  the 
square  of  the  error.  The  objective  can  be  stated  as 

Minimize  =  (2.3) 

m=l 
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where,  is  the  voltage  at  the  electrode  predicted  by  the  model  while  is  the 
voltage  measured  at  the  same  electrode  experimentally  and  is  the  number  of  data 
points  or  electrodes  at  which  voltages  are  measured.  The  variables  of  the  optimization 
problem  are  the  resistivity  values  in  the  three  principal  directions  that  are  the 
components  of  the  vector  p .  The  resistivity  values  are  used  as  variables  instead  of  the 
conductivity  values  because  the  voltages  have  an  inversely  proportional  relationship 
with  the  conductivities  that  would  make  the  objective  function  highly  nonlinear. 

The  Gauss-Newton  algorithm  linearly  approximates  the  error  at  each  iteration  and 
then  minimizes  the  resultant  problem  using  the  Newton's  method.  Therefore,  this 
approach  requires  only  the  computation  of  the  gradient  of  the  error.  The  application  of 
this  approach  to  compute  resistivity  by  minimizing  the  error  in  the  computed  voltages 
is  summarized  below.  The  gradient  of  the  objective  function  defined  in  equation  (2.3), 
can  be  expressed  as 


dF 

Qpi 


nd  ,  . 

m=\ 


Ml 

dPi 


(2.4) 


Using  Newton's  method,  the  optimality  criterion,  VF  =  0  can  be  solved  iteratively 
by  updating  the  resistivity  values  at  the  fc*  iteration  as 


^  k+\  ,  \  ^k 

Pi  ^Pi  +^Pi 


(2.5) 


Where  the  update  vector  Ap/  is  computed  by  solving 


dp,  %  dp, dp, 


(2.6) 


The  second  derivative  of  the  objective  function  or  the  Hessian  matrix  can  be 
computed  as 


Hij  =  ^  ^ 

dp,dpj 


=  21 


54  54 


+ 


=ii  Spi  Spj  SpiSpj 


(4-<C) 


(2.7) 


The  Gauss-Newton  approach  involves  approximating  the  error  (or  in  this  case 
voltage)  as  linear  at  each  iteration  so  that  the  second  derivative  of  the  voltage  is  set  to 
zero,  approximating  the  Hessian  matrix  as 


=2N 

;ti[^pi'spjj 


(2.8) 
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Using  this  approximation,  equation  (2.6)  is  solved  to  compute  the  resistivity 
update  vector  .  In  order  to  implement  this  iterative  strategy  to  compute  the 
resistivity  values,  it  is  necessary  to  compute  the  gradient  of  the  voltage  with  respect  to 
the  resistivity.  This  gradient  is  needed  for  computing  both  the  gradient  and  the  Hessian 
of  the  objective  function.  It  can  be  computed  by  taking  the  gradient  of  the  governing 
equations.  The  finite  element  method  converts  the  governing  equation  (2.1),  into  a  set  of 
linear  simultaneous  equations,  often  express  in  the  form. 

Where  [K]  is  the  global  conductivity  matrix,  {O}  is  the  voltage  vector  containing 
the  nodal  values  of  the  potential  or  voltage  and  the  current  vector  {1}  contains  the 
contribution  from  current  applied  at  the  boundaries.  Taking  the  derivatives  of  both 
sides  of  this  equation  with  respect  to  the  variables  ,  we  get. 


{®} 


(2.10) 


The  current  sources  are  clearly  not  a  function  of  the  resistivity  values  and 
therefore  the  current  vector's  derivative  is  zero.  The  right  hand  side  of  equation  (2.10) 
can  be  computed  element  by  element  and  assembled  to  create  a  global  vector.  Then  the 
gradient  of  all  the  nodal  voltages  with  respect  to  the  resistivity  values  can  be  obtained 
by  solving  equation  (2.10).  Of  course,  we  need  the  gradient  only  for  the  nodes  that  are 
located  at  the  electrodes  where  the  voltages  are  measured. 


5.  Numerical  validation 

The  inverse  approach  described  in  the  previous  section  was  implemented  on  a 
finite  element  program  for  solving  electrostatic  problems.  In  order  to  first  verify  the 
validity  of  the  algorithms  and  the  implementation,  the  model  itself  was  used  to  create 
data  by  computing  the  voltages  at  the  electrodes  using  assumed  values  of  resistivity. 
Then  using  this  data  in  the  inverse  approach,  the  resistivity  was  computed  starting  from 
random  values  to  verify  if  the  known  correct  value  can  be  computed.  A  plate  whose 
dimensions  are  54x52x1.5  was  modeled  for  the  validation.  The  total  number  of 
elements  in  the  TEA  model  was  3024  where  each  element  is  of  size  2x2x0.5  mn^ .  A 
constant  current  of  28.9mA  was  applied  at  the  electrodes.  For  generating  the  data  used 
for  inverse  algorithm  validation,  the  material  was  modeled  as  isotropic  with  a 
resistivity  of  0.01  mClm  in  all  three  principal  directions.  Using  the  voltage  data 
computed  using  this  model,  the  inverse  problem  was  solved  to  estimate  resistivity 
values  and  in  this  case  exact  solution  were  obtained  for  the  resistivity.  Therefore, 
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random  noise  from  0.1  to  1  per  cent  was  added  to  the  voltage  data  to  simulate  typical 
experimental  error  in  measured  data.  The  inverse  problem  was  solved  repeatedly  for 
data  with  different  noise  levels  and  the  effect  of  the  added  noise  was  evaluated  to 
determine  the  resultant  error  in  the  resistivity  computed  by  the  inverse  approach. 

Fig.  5  shows  a  plot  of  the  error  in  the  computed  resistivity  values  versus  the 
percentage  random  noise  or  error  introduced  into  the  voltage  data.  For  data  with  no 
error,  all  three  principal  resistivity  values  are  computed  with  zero  error  but  as  the  noise 
level  increases  the  errors  in  the  computed  values  increase.  With  0.1%  random  noise 
added  to  the  voltage  data,  the  errors  in  the  resistivity  values  range  from  0.03%  to  0.17%. 
In  the  case  of  0.5%  random  noise,  the  errors  are  from  0.38%  to  0.73%  and  for  1%  noise, 
they  are  from  0.85%  to  5.39%.  Clearly,  the  computed  values  of  the  resistivity  are 
sensitive  to  errors  in  the  data.  The  resistivity  in  the  thickness  was  found  to  be  the  most 
sensitive  to  the  noise. 


Fig.  5  Error  in  computed  resistivity  due  to  noise  in  voltage  data 

It  the  material  does  not  have  uniform  properties,  then  an  average  value  is  obtained 
using  the  inverse  approach.  Therefore,  if  the  material  has  a  region  with  different 
properties  due  to  embedded  inclusions  or  due  to  damage  /  delamination  then  the 
computed  average  value  would  be  significantly  different  than  the  undamaged  uniform 
material's  properties.  This  can  serve  as  an  indicator  for  detecting  defects  or  monitoring 
damage  in  composite  structures.  To  simulate  this,  models  with  varying  size  regions 
with  different  properties  were  modeled.  Fig.  6  shows  models  with  square  regions  that 
represent  damaged  regions  that  have  different  properties  than  the  surrounding 
material. 
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(D)  22  X  22 


(A)4x4  (B)  10x10  (C)  16x16 

Fig.  6  Size  and  shape  of  the  damaged  region  modeied  in  the  forward  problem 

All  the  models  are  plates  of  dimensions  50x50x1  mm^,  which  are  assumed  to  be 
made  of  woven  composites  whose  conductivity  values  are  o-j  =(T2  =15  S/mm  and 
<73  =  1.5  S/mm.  The  damaged  region  is  assumed  to  have  very  low  conductivity  values: 
cTj  =  0-2  =  (7,  =  ixi0“‘  S/mm.  Using  these  models,  voltage  data  at  the  electrodes  where 
generated  to  be  used  as  data  for  testing  using  the  inverse  approach.  Fig.  7  shows  the 
average  values  of  conductivity  computed  by  the  inverse  method.  The  largest  change  in 
conductivity  occurs  in  the  through  thickness  direction  in  this  case  where  we  have 
assumed  that  the  damaged  region  is  isotropic  with  very  low  conductivity. 


Fig.  7  Computed  average  conductivity  versus  damage  size 
6.  Experimental  Results 

Several  plate-like  specimens  made  of  both  unidirectional  and  woven  composites 
were  used  to  make  experimental  measurements  of  the  voltages  on  electrodes  located  on 
both  sides  of  the  plate.  The  voltage  data  measured  from  experimental  specimen  was 
used  to  determine  the  resistivity  values,  first  for  undamaged  composite  plate  like 
specimen  with  uniform  values  of  resistivity.  Thereafter,  specimen  with  artificial  damage 
or  non-uniform  properties  was  used  to  determine  an  average  value  of  resistivity. 
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Resistivity  measurement  for  undamaged  unidirectionai  composite  specimen 

Continuous  carbon  fiber  unidirectional  composite,  with  35%  resin  content  and 
eleven  layers,  was  used  to  prepare  specimens  of  two  different  size,  which  are 
52x47x1.5  for  the  first,  and  54x52x1.5  for  the  second.  Both  specimen  had  eight 
electrodes,  each  about  1x1  mrn^ ;  four  electrodes  on  the  top  surface  and  the  other  four 
electrodes  on  the  bottom  surface.  The  unidirectional  composite  specimens  were  made 
by  cutting  it  out  of  commercially  available  composite  sheets.  In  order  to  make  the 
electrodes,  the  electrode  area  was  first  polished  with  sandpaper  to  remove  the  surface 
layer.  Thereafter,  silver  paint  was  applied  to  this  area  and  after  it  was  dried,  a  thin 
copper  strip  was  attached  to  this  areas  using  silver  epoxy  to  adhere  to  copper  strip  and 
to  provide  good  electrical  contact.  The  measurement  probes  where  attached  to  the 
copper  strips  for  measuring  voltages  and  for  injecting  currents.  The  eight  electrodes 
were  numbered  as  shown  in  Fig.  3  . 

Table  I:  Measurement  configurations 


Current  Pairs  Voltage  measurement  pairs  (Target,  Reference) 

(Source,  Drain) _ 


(3,1) 

(2,1) 

(4,1) 

(5,1) 

(6,1) 

(7,1) 

(8,1) 

(4,1) 

(2,1) 

(3,1) 

(5,1) 

(6,1) 

(7,1) 

(8,1) 

(7,1) 

(2,1) 

(3,1) 

(4,1) 

(5,1) 

(6,1) 

(8,1) 

(8,1) 

(2,1) 

(3,1) 

(4,1) 

(5,1) 

(6,1) 

(7,1) 

In  table  I,  the  electrode  pairs  at  which  current  was  injected  are  listed  and  for  each 
such  pair  voltage  was  measure  at  six  different  electrode  pairs.  Voltage  difference  is 
measured  by  a  Keithley  2002  multimeter  using  the  four-probe  method  and  DC  constant 
current  magnitude  of  7.235mA  was  injected  at  the  current  pairs. 

Uncertainty  of  voltage  measurement  is  about  0.0026%  using  Keithley  2002.  In 
order  to  calculate  the  measurement  error,  20  consequent  voltage  data  were  taken  on 
each  voltage  measurement  pair  to  calculate  the  average  value  and  the  standard 
deviation.  Numerical  validation  results  suggest  that  if  the  error  in  voltage  data  is 
0.0026%,  the  calculated  resistivity  values  should  have  less  0.1%  error.  However,  in 
practice,  the  error  can  be  higher  due  to  inaccuracy  in  modeling  the  geometry  of 
specimen,  electrodes  and  also  numerical  errors  in  solving  the  forward  problem.  The 
calculated  values  of  resistivity  are  an  average  value  for  the  entire  specimen,  which  is 
assumed  to  have  uniform  properties  in  the  model.  Two  different  geometrically  accurate 
FEA  models  were  created  corresponding  to  the  two  specimens  used  in  the  experiments. 
The  FEA  model  for  the  first  specimen  has  10176  elements  with  the  element  size  of 
lxlxO.5  .  The  second  specimen  has  11660  with  the  same  element  size. 

Table  II  shows  computed  values  of  resistivity  in  each  direction  for  the  two 
specimens.  For  both  specimen,  the  electrodes  where  modeled  as  current  injection  areas 
of  1x1  mm^ .  In  the  experimental  specimen,  it  is  difficult  to  create  accurately  sized 
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electrodes.  To  study  the  sensitivity  of  the  results  to  error  in  modeling  the  electrode  area, 
the  same  voltage  data  was  used  with  models  that  used  2x2  mrn^  current  injection  regions 
to  model  the  electrodes.  The  computed  values  of  resistivity  did  not  change  significant 
indicating  that  it  is  not  necessary  to  model  the  electrodes  with  great  precision.  The 
resistivity  values,  in  table  II,  show  that  the  resistivity  in  the  transverse  and  thickness 
direction  are  several  orders  of  magnitude  larger  that  the  resistivity  in  the  fiber  direction. 


Table  II:  Computed  values  of  resistivity 


Direction 

52x47x1.5  specimen 

54x52x1.5  specimen 

Fiber  direction 

0.024  mQm 

0.020  mQm 

Transverse  direction 

18.6  mQ.m 

13.4  mOm 

Thickness  direction 

67.5  mClm 

60.0  mClm 

Least  square  error: 

5.25x1 0'^^ 

3.68x10'^ 

Woven  composite  specimen  with  and  without  damage 


Resistivity  values  were  determined  for  woven  composite  specimen  also  using  the 
same  procedure.  To  study  the  effect  of  damage,  artificial  damage  was  introduced  by 
embedding  Teflon  patches  between  the  layers.  We  fabricated  both  the  damaged  and 
undamaged  specimen  using  four  layers  of  woven  carbon  fiber  prepregs.  To  create  the 
damaged  specimen.  Teflon  patches  were  introduced  between  the  first  and  second  layers 
as  well  as  the  third  and  fourth  layers  as  shown  in  the  Fig.  8  .  The  specimen  size  was 
50x50x1  mm^  and  the  size  of  the  Teflon  patches  were  10x10  mm^  and  they  were 
centrally  located  within  the  specimen.  Eight  electrodes  attached  at  the  corners  of  the 
specimen  were  used  to  gather  voltage  data.  The  electrodes  were  created  by  inserting 
copper  strips  between  the  first  and  second  layers  as  well  as  between  the  third  and 
fourth  layers  of  prepreg  before  curing.  This  method  avoids  the  need  for  attaching 
electrodes  with  silver  epoxy  and  it  provides  better  contact  as  well  as  robust  connection 
that  are  not  easily  damaged. 


layer 
3*^**  layer 


Embedded  Teflon  Patch 


(a)  Section  showing 
layers 


Fig.  8  Woven  fiber  composite  prepregs  with  embedded  Tefion  patches 
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To  determine  the  resistivity,  three  undamaged  specimen  and  three  Teflon 
embedded  specimen  were  used.  The  results  obtained  are  listed  in  Table  III,  which 
shows  values  for  resistivity  for  the  six  specimen  and  the  average  values.  Clearly,  there 
is  significant  change  in  resistivity  due  to  the  Teflon  embedding  that  can  be  detected 
using  the  resistivity  measurement  technique  described  here.  This  approach  therefore 
has  potential  to  be  used  for  quality  control,  to  detect  manufacturing  defects  such  as 
voids  and  air  gaps  as  well  as  delamination  or  damage  that  may  occur  during  usage. 


Table  III:  Resistivity  values  for  woven  fiber  composite  specimen 


Specimen  No. 

Specimen  without  teflon 

Specimen  with  Teflon 

1 

CT^  =  15.39[S/mm] 

o-;  =  9.4086[S/mm] 

=  17.25[S/mm] 

cr'  =11. 88[S/mm] 

0-3  =  1.61 7[S/mm] 

0-'  =  0.9708[S/imn] 

2 

CTj  =  13.46[S/mm] 

o-;  =  9.2415[S/imn] 

G^  =  18. 15  [S/mm] 

cr'  =  8.539[S/mm] 

0-3  =1.18[S/mm] 

0-'  =  0.7438[S/imn] 

3 

G^  =  22.86[S/mm] 

a-;  =  8.2043[S/mm] 

Gr^  -  22.32 [S/mm] 

a-'  =  8.9204[S/mm] 

0-3  =  1.02[S/mm] 

cr'  =  0.8232[S/mm] 

Average 

o-j  =  17.24[S/imn] 

0-;  =8.9515[S/mm] 

0-2  =  19.24[S/mm] 

cr'  =  9.7798[S/mm] 

CTj  =  1.272[S/mm] 

cr'  =  0.8459[S/mm] 

7.  Discussion 

An  algorithm  for  solving  inverse  problem  to  determine  average  resistivity  values 
in  composite  structures  was  demonstrated.  The  method  can  use  data  from  an  arbitrarily 
large  number  of  electrodes  to  compute  average  values  of  resistivity  or  conductivity  for 
the  structure.  Finite  element  models  for  the  structure  are  used  to  solve  the  forward 
problem,  making  this  method  very  general  and  applicable  to  arbitrary  shaped 
structures.  Ideally,  the  electrodes  should  be  embedded  in  the  structures  during  the 
manufacturing  process  itself  so  that  it  can  be  used  for  quality  control,  detection  of 
defects  as  well  as  subsequent  health  monitoring.  One  of  the  advantages  of  measuring 
resistivity  is  that  damage  can  be  detected  even  in  structures  that  were  not  tested  during 
manufacturing.  Damage  can  be  detected  for  structures  that  are  in  use  by  attaching 
electrodes  on  the  surface,  determining  the  average  resistivity  and  comparing  it  to  values 
associated  with  undamaged  material.  The  main  source  of  error  in  this  approach  arises 
from  inaccuracy  in  the  geometric  models  of  the  structure  and  the  electrodes.  Random 
noise  added  to  the  voltage  data  used  in  numerical  validation  indicates  that  any  error  in 
the  voltage  data  can  get  amplified  due  to  difficulties  in  numerical  convergence. 
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In  principle,  a  similar  inverse  approach  can  also  be  used  to  determine  applied 
loads  on  the  structure.  Preliminary  experimental  studies  indicate  that  this  may  be 
difficult  because  the  changes  in  the  electrode  voltage  due  to  strains  can  be  very  small. 
Even  with  amplification,  the  data  is  hard  to  use  because  of  significant  non-linearity  in 
the  observed  behavior.  However,  the  approach  we  developed  is  promising  for  detecting 
damages  or  defects  because  they  cause  significant  changes  in  resistivity  and  is  therefore 
easier  to  detect.  Further  study  is  needed  to  explore  ways  of  determining  applied  loads 
and  strains. 
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III.  COMPOSITE  ACTUATING  STRUCTURES 


1.  Overview 

Polymer  matrix  composites  are  widely  used  as  a  structural  material  in  a  variety  of 
aerospace  applications  including  commercials  jets,  unmanned  air  vehicles  as  well  as 
micro  air  vehicles.  In  many  of  these  applications  weight  is  an  important  constraint  and 
available  space  is  very  limited.  This  provides  an  incentive  to  pack  as  much  functionality 
as  possible  into  the  structure  itself  by  making  it  multifunctional.  For  air  vehicle 
applications,  structures  with  actuation  ability  are  particularly  desirable  for  designing 
active  structures  such  as  morphing  wings/body  panels  or  flapping  wing  like  structures. 
Methods  of  actuation  studied  in  the  past  have  mainly  consisted  of  embedded  peizo- 
electric  fibers  and  patches.  In  this  project,  electromagnetic  means  of  actuation  were 
studied  wherein  embedded  ferromagnetic  materials  are  used  for  the  actuation  using 
external  or  internally  generated  magnetic  fields. 

Composite  structures  made  of  epoxy  matrix  and  carbon  fiber  reinforcement  have 
excellent  structural  properties  including  stiffness/rigidity  and  high  strength.  Therefore, 
many  aircraft  structures  are  made  of  such  composites,  particularly  in  unmanned  and 
micro  air  vehicles.  Due  to  the  conductivity  of  carbon  fibers,  it  is  possible  to  conduct 
currents  through  these  composites.  Magnetic  field  can  be  generated  in  composite 
structures  by  current  flowing  though  the  reinforcing  fiber  s/conducting  wires  as  well  as 
due  to  embedded  permanent  magnets.  The  magnetic  forces  generated  on  the  structure 
can  be  large  enough  to  cause  deformation  if  ferromagnetic  materials  are  embedded  in 
the  structure  and  large  currents  are  flowing  through  the  structure.  Unlike  in  traditional 
electrical  machinery,  structures  that  are  meant  to  actuate  are  designed  to  deform 
significantly  due  to  the  magnetic  forces.  Therefore  the  structure  cannot  be  treated  as  a 
rigid  body  and  its  deformation  needs  to  be  computed  using  an  analysis  model  that 
couples  the  magnetic  and  structural  models.  To  compute  the  resultant  deformations, 
strains,  and  stresses,  a  two  stage  analysis  approach  is  adopted  here.  The  magnetostatic 
problem  is  solved  first  to  compute  the  magnetic  flux  density  and  field  distribution.  This 
result  is  then  used  to  compute  the  body  forces  generated  on  the  structure  due  to 
magnetic  forces.  In  the  second  stage  a  solid  mechanics  analysis  is  performed  using  these 
magnetic  forces  to  compute  the  mechanical  deflections.  In  general  if  these  deflections 
are  large,  causing  the  current  carriers  and  embedded  soft  and  hard  magnets  to  move 
significantly  relative  to  each  other,  then  the  movement  would  alter  the  magnetic  field, 
requiring  an  iterative  solution.  However,  in  this  project,  only  linear  elastic  deformation 
with  a  two  stage  analysis  as  described  was  considered.  For  the  type  of  actuators  that  is 
most  appropriate  for  this  application,  linear  models  are  sufficient  as  described  in  later 
chapters.  The  structure  is  also  modeled  as  linear  elastic  since  the  deformation  of  the 
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mechanism  is  reversible.  It  would  be  beneficial  in  future  work  to  incorporate  large 
deformation  models  to  simulate  the  full  range  of  motion  of  highly  flexible  structural 
mechanisms. 


2.  Motivation  for  structurai  actuators 

The  primary  application  that  serves  as  motivation  for  this  work  is  the  design  of 
micro  air  vehicles  (MAV)  and  unmanned  air  vehicles  (UAV).  The  goal  was  to  develop 
computational  tools  for  design  and  analysis  of  magnetically  actuated  structural 
mechanisms  for  these  air  vehicles  so  that  structural  components  of  the  vehicle  including 
wings,  fuselage  and  tail  can  be  designed  function  as  structural  actuators.  Therefore, 
rather  than  have  external  mechanisms,  linkages,  motors  or  actuators  to  produce  the 
necessary  motions,  the  structure  is  designed  to  deform  in  specified  manner  due  to  built- 
in  magnetic  actuation  capabilities.  This  would  allow  these  structures  to  flex  or  change 
shape  for  morphing  applications  as  well  as  oscillate  or  vibrate  to  produce  flapping 
motion. 

Fig.  9  shows  the  conceptual  design  of  a  simple  flapping  wing  actuation 
mechanism  where  the  wings  are  attached  to  a  flexible  structural  support  that  is  actuated 
by  the  electromagnets  built-into  this  support.  The  key  design  challenges  for  this  concept 
include: 

(i)  Designing  magnetic  actuators  that  are  strong  enough  to  produce  the 
necessary  deformation 

(ii)  Designing  the  support  structure  and  wing  such  that  the  desired  motion  is 
produced  when  the  actuating  magnets  are  activated 

(iii)  Designing  the  dynamics  of  the  wing  so  that  at  resonance  the  flapping 
motion  will  produce  the  desired  mode  of  vibration  that  can  generate  thrust 
and  lift. 
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Fig.  9  Flapping  wing  conceptual  design 

The  actuator  can  be  better  integrated  with  the  structure  if  the  deforming 
component  is  built  as  part  of  the  wing  structure.  A  design  based  on  this  concept  is 
shown  in  Fig.  10  . 


Fig.  10  Alternate  design  of  actuator 

The  structure  is  designed  to  be  compliant  so  that  it  can  deform  to  produce  the 
desired  actuation.  The  idea  is  illustrated  using  a  simple  wing  design  in  Fig.  11  where 
the  magnetic  forces  are  shown  as  acting  along  the  edge  of  half  of  the  wing  assembly. 
The  force  causes  the  wing  structure  to  deform  to  produce  the  flapping  motion.  The 
actual  deformation  mode  will  depend  not  only  on  the  shape  of  the  wing  and  the 
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support  but  also  on  the  reinforcements,  fiber  direction  and  number  layers  in  the 
composite  shell. 


Fig.  11  Deforming  structurai  mechanism 

In  order  to  design  the  embedded  actuators,  it  is  necessary  to  have  the  ability  to 
compute  the  magnetic  field  produced  by  the  actuator,  compute  forces  on  structures  and 
then  compute  the  deformed  shape  of  the  structure  due  to  these  forces.  To  design  the 
structure  itself,  we  would  like  to  compute  the  shape  of  the  support  structure,  the 
orientation  of  the  fibers  of  the  composite,  the  possible  location  of  holes  or 
reinforcements  such  that  the  structure  would  deform  or  oscillate  in  the  desired  mode.  A 
possible  design  tool  for  computing  the  geometry  and  reinforcement  is  topology 
optimization.  This  idea  was  explored  as  part  of  this  project  and  some  of  the  results  are 
presented  in  the  next  section. 

3.  Designing  the  shape  and  topology  of  structurai  actuator 

In  order  to  design  the  shape  of  the  structure  that  is  appropriate  for  the  desired 
actuation,  the  geometry  design  problem  is  stated  as  a  design  optimization  problem. 
Firstly,  a  region  within  which  the  geometry  must  fit  is  defined  as  a  feasible  region.  The 
geometry  is  defined  within  this  region  as  the  level  set  or  contour  of  a  density  function 
(or  the  characteristic  function).  Contours  of  this  function  corresponding  to  a  threshold 
value  are  defined  as  the  boundaries  of  the  shape  so  that  regions  where  the  value  of  the 
function  is  below  the  threshold  are  not  part  of  the  geometry.  Hence,  the  boundary  may 
be  defined  using  the  following  implicit  equation 


<t>(x,y)-(t),h  =0 


(3.1) 
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Fig.  12  Shape  representation  using  shape  density  function 

In  traditional  topology  optimization  methods,  porosity  or  density  of  the  material 
is  treated  as  constant  within  each  element.  In  this  work,  the  density  was  assumed  to 
vary  continuously  within  the  feasible  region.  Fig.  12  illustrates  this  shape  representation 
where  the  rectangle  represents  the  feasible  region,  the  arrows  at  the  top  represent  a 
uniformly  distributed  load  and  the  structure  is  to  be  supported  at  the  bottom.  A  mesh  is 
generated  for  the  feasible  region  and  the  density  function  is  defined  within  this  feasible 
region  by  piece-wise  interpolation  within  elements  of  the  mesh.  Contours  of  the  density 
function  are  plotted  in  the  figure.  The  contour  of  the  density  function  corresponding  to 
the  threshold  value  (|)^  is  the  boundary  of  the  solid  and  the  regions  with  higher  values 
of  density  is  the  interior  of  the  solid  shown  as  the  shaded  region  in  the  Fig.  12  .  In  this 
example,  there  are  multiple  contours  that  correspond  to  the  threshold  value  each 
representing  part  of  the  boundary.  Shape  representation  using  a  contour  of  the  density 
function  enables  the  entire  geometry  to  be  treated  as  a  variable.  By  changing  the  density 
function  it  is  possible  to  not  only  modify  existing  boundaries  but  also  to  create  new 
internal  boundaries. 

The  density  function  (|)  =  1 ,  where  the  material  is  fully  dense  and  (|)  <  (|),^  where  there 
is  no  material.  The  density  function  can  also  take  on  intermediate  values  but  as 
explained  later  the  relation  between  density  and  material  properties  is  selected  such 
that  the  optimal  designs  are  close  to  fully  dense.  A  new  internal  boundary 
corresponding  to  a  hole,  for  example,  would  be  created  if  the  value  of  the  density 
function  decreases  to  the  threshold  value  in  a  region. 

To  define  the  density  function,  the  feasible  domain  is  divided  into  triangular  or 
rectangular  elements.  The  density  function  is  interpolated  within  each  element.  A 
contour  of  the  density  function  corresponding  to  the  threshold  value  passes  through  the 
element  if  some  nodes  of  the  element  have  nodal  density  values  higher  than  the 
threshold  value  while  others  have  values  below.  In  this  project,  linear  or  bilinear 
elements  as  well  as  B-spline  elements  were  used.  The  contour  is  therefore  plotted  by 
joining  the  points  along  the  contour  that  has  density  value  equal  to  the  threshold  value. 
A  C°  continuous  density  function  ensures  C°  continuous  boundaries  for  the  final  shape. 
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whereas,  quadratic  and  cubic  B-spline  elements  provide  a  and  continuous  density 
functions  and  contours  respectively. 

The  mesh  used  for  defining  the  density  function  can  also  serve  as  the  finite 
elements  for  sfructural  analysis.  The  values  of  fhe  density  function  at  the  nodes  serve  as 
the  design  variables  of  the  optimization  problem.  Initially,  the  nodal  density  values  are 
set  equal  to  unity  for  all  the  nodes  so  that  the  geometry  is  identical  to  the  feasible 
region.  During  fhe  optimization  process  the  nodal  density  values,  (|)i  are  modified  by 
the  optimization  algorithm  which  iteratively  searches  for  the  optimal  values  of  the 
nodal  densities  such  that  the  defined  objective  function  is  minimized. 

For  the  application  of  interest  in  this  project,  the  objective  is  to  design  a  structure 
that  deforms  in  a  particular  fashion  when  subjected  to  electromagnetic  forces.  Therefore 
the  objective  function  is  defined  as  the  error  between  the  desired  deflection  and  the 
actual  deflection.  In  this  project,  the  method  was  applied  only  to  planar  (2D)  problems, 
such  as  plane  stress  and  plane  strain.  In  addition  to  minimizing  this  objective  function, 
a  constraint  on  the  total  mass  of  the  structure  is  applied. 

Structural  synthesis  is  the  inverse  of  the  structural  analysis  problem.  The 
structural  analysis  problem  is  typically  stated  as  a  principle  of  virtual  work  (PVW), 

j  {8e}‘  [D]  {e}  dfi  ==  •  8u  dn  +  j  •  8u  dr  (3.2) 

Q  D  Tj 

The  domain  Q  represents  the  shape  and  topology  of  the  component  whose 
structural  properties  are  being  analyzed.  The  finite  element  method  is  used  to  solve  for 
the  displacement  field  u(x)  for  every  x  belonging  to  Q. 

The  shape  and  topology  synthesis  problem  involves  solving  for  a  domain  Q  thaf 
optimizes  some  structural  property  for  given  loading  and  boundary  conditions.  The 
geometry  Q  is  defined  as  the  region  within  the  original  feasible  domain  Qq  where  the 
shape  density  function  has  a  value  greater  than  the  threshold  value.  The  minimization 
problem  can  be  stated  as 


m 

Minimize  !!((())  =  ^(U;  -  u* 

i=0 


(3.3) 


subject  to. 

M(^)  =  1  ^  dQ  <  Mo 

Q. 

(3.4) 

j  {58}‘  [Dm  {£}  dQ 

(3.5) 
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(3.6) 


n((|))  is  the  sum  of  the  square  of  errors  in  the  computed  nodal  displacement  at 
nodes  where  a  desired  value  of  displacement  has  been  specified,  where,  (|)(x)  is  the 
density  function.  It  is  assumed  that  a  desired  value  of  displacement  u*  has  been 
specified  at 'm'  nodes  and  U;  is  the  computed  value  of  the  displacement.  Equation  (3.4) 
describes  the  constraint  that  the  mass  M  of  the  component  should  be  less  than  or  equal 
to  a  specified  value  Mo. 

The  optimization  problem  stated  in  equations  (3.3)-(3.6)  can  be  solved  using 
mathematical  programming  techniques  or  optimality  criteria  methods.  The  objective 
function  is  non-linear  and  the  constraint  on  weight  is  linear.  Each  evaluation  of  the 
objective  function  requires  a  computationally  expensive  finite  element  analysis  to 
compute  the  displacement  at  the  nodes  of  the  finite  element  mesh  for  the  structure. 
Therefore,  an  algorithm  that  does  not  require  excessive  number  of  function  evaluations 
is  preferred.  A  modified  form  of  sequential  linear  programming  [29]  was  used  for  the 
results  presented  in  this  report. 

When  the  shape  defined  by  the  density  function  varies,  the  structural  properties 
must  vary  accordingly.  This  implies  that  the  material  property  coefficients  defined  in 
the  matrix  [D]  must  depend  on  the  density  function  (|)(x,y).  We  seek  relations  that  are 
simple  and  therefore  easy  to  integrate  over  each  element  when  density  varies  linearly 
within  each  element.  In  addition  we  would  like  relations  that  lead  to  clearly  defined 
topologies  so  that  the  final  shape  obtained  is  fully  dense  and  the  density  function 
transitions  sharply  at  the  boundary  from  full  density  to  the  lowest  possible  (threshold 
value). 

The  material  property-density  relation  should  be  such  that  if  the  density  decreased 
in  a  region,  the  stiffness  should  decrease  causing  the  material  to  become  weaker  in  that 
region.  This  would  be  achieved  if  the  slope  of  the  objective  function  with  respect  to  the 

variables,  -^n(u),  is  negative.  The  optimization  algorithm  would  therefore  decrease  the 

4S,. 

density  in  regions  where  material  is  under-utilized  causing  either  new  boundaries 
(holes)  to  be  created  in  such  regions  or  causing  existing  boundaries  to  shrink  inwards. 
We  have  used  some  linear  and  non-linear  material  property  density  relations  that 
satisfy  this  criterion. 

Homogenization  method  has  been  used  to  determine  the  relation  between 
porosity  (or  density)  and  elastic  constants  by  assuming  a  microstructure.  Typically  a 
square  unit  cell  with  a  circular,  square  or  rectangular  void  is  used  to  determine  the 
elastic  constants.  The  size  of  the  void  is  changed  to  vary  the  porosity  (or  density)  and 
the  elastic  constant  are  computed  for  various  values  of  porosity.  Since  the 
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homogenization  process  is  computationally  expensive  the  elastic  constants  are 
computed  for  a  few  values  of  porosity  and  then  a  curve  is  fitted  over  these  computed 
points  to  obtain  a  relation.  However,  the  relation  obtained  is  not  unique  because  it 
depends  on  the  microstructure  assumed.  The  optimal  designs  that  are  obtained  will  also 
be  therefore  different  based  on  the  assumption  used.  This  raises  the  question  as  to 
which  relation  is  ideal  for  computing  optimal  shapes.  Since  the  real  material  does  not 
have  varying  density  or  porosity  it  is  preferred  that  the  optimal  designs  are  fully  dense. 
Any  relation  that  leads  to  such  design  is  therefore  preferable. 

In  our  implementation,  polynomial  relation  was  assumed  between  elastic  constant 
and  the  density  function.  For  example,  it  can  be  assumed  that  the  elastic  modulus  of  the 
material  is  a  quadratic  function  of  the  density.  Similarly,  one  could  use  higher  order 
approximations.  Just  as  different  microstructure  assumptions  lead  to  different  optimal 
designs  when  homogenization  method  is  used,  different  polynomial  relations  between 
elastic  modulus  and  density  lead  to  different  designs.  The  criterion  that  we  used  for 
selecting  the  relation  is  the  sharpness  of  the  density  transition  at  the  boundary.  In  other 
words,  we  want  the  material  inside  the  shape  to  be  fully  dense  ((|)=1)  and  the  material  to 
have  the  lowest  possible  density  where  the  holes  are  located.  At  the  boundary  we  want 
the  density  to  transition  sharply  from  the  highest  value  to  the  lowest  value,  so  that  we 
have  clear  and  well  defined  boundaries.  When  a  linear  relation  is  assumed  sharp 
boundaries  are  not  obtained  except  when  the  threshold  value  is  set  close  to  1.  It  was 
found  that  in  general,  higher  order  approximations  of  the  material  property-density 
relations  lead  to  the  desired  behavior. 

Assuming  order  polynomial  relation  between  the  elasticity  modulus  and  the 
density  function,  we  get  the  following  material  property-density  variation  for  plane 
stress  problems. 


d„  = 


d,2  = 


Ef 

1-v" 

Eyf 

1-v" 

Ef 

2(1 +  v) 


(3.7) 


The  coefficients  djj  are  the  elements  of  the  elasticity  matrix.  Note  that  in  the  above 
relation  we  do  not  assume  that  Poisson's  ratio  changes  with  density  and  therefore,  for 
this  approximation  the  material  coefficients  reduce  to  zero  as  density  goes  to  zero,  that 
is,  dij=0  for  (|)=0.  The  elasticity  matrix  can  then  be  conveniently  defined  as 

[Dp(^)]  =  [D]f  (3.8) 
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where,  [D]  is  the  elasticity  matrix  for  plane  stress  or  plane  strain. 

In  order  to  use  a  mathematical  programming  algorithm  to  compute  the  optimal 
design,  it  is  necessary  to  compute  the  gradient  of  the  objective  function  and  the 
constraint.  The  gradient  of  the  objective  function  is: 


an 


'Vor  ».aUj 
5(|)i 


(3.9) 


The  gradient  of  nodal  displacements  can  be  computed  using  the  standard  design 
sensitivity  analysis  methods.  The  equilibrium  equations  are  reduced  to  a  set  of  linear 
simultaneous  equations  by  the  finite  element  method  which  is  usually  expressed  as 

[K]{U}  =  {F}  (3.10) 


{U}  is  the  displacement  vector  that  contains  Uj,  the  displacement  components  at 
the  nodes  and  {F}  is  the  load  vector.  The  gradient  of  Uj  can  be  computed  by  solving  the 
equation 

[K]^  =  W  (3.1.) 

where,  {?i}  is  the  adjoint  variable  that  is  computed  as 

W  =  -^(U}  (3.12) 


4.  Examples  of  Topology  Design 
Gripper  Mechanism  Design 

As  an  example,  let  us  consider  the  design  of  a  mechanical  gripper.  The  mechanism 
is  supported  at  the  two  corners  on  its  left  edge  and  input  forces  of  magnitude  50000  N 
are  applied  in  the  middle  of  the  left  edge  as  shown  in  Fig.  13  .  Vertical  displacements 
desired  at  the  points  A  and  B  that  causes  them  to  move  towards  each  other  or  to  grip  a 
work-piece.  Forces  of  magnitude  5000  N  are  applied  at  the  points  A  and  B,  where 
displacements  are  expected,  to  model  the  resistance  of  the  work-piece  once  the 
mechanism  comes  in  contact  with  the  work-piece. 

The  size  of  the  design  domain  is  5  x  5  m  as  shown  in  Fig.  13  .  Displacements  of 
magnitude  0.00025  are  prescribed  at  two  corner  points  A  and  B  and  displacements  of 
magnitude  0.000025  are  prescribed  at  the  points  where  the  input  forces  are  applied. 
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Fig.  13  Feasible  domain  for  gripper  mechanism  design  with  a  30  x  30  mesh 

The  material  of  the  domain  is  assumed  to  be  steel  with  modulus  of  elasticity  equal 
to  200  GPa  and  the  Poisson's  ratio  of  0.3.  The  original  domain  has  been  discretized  with 
a  sparse  mesh  of  30  x  30  elements. 


(a)  (b)  (c) 


Fig.  14  Topology  results  for  a  mechanical  gripper  design  with  a  30  x  30  mesh  using  (a) 
Quad  4N  elements  (b)  B-spline  9N  elements  (c)  B-spline  16N  elements 

The  topology  results  of  the  optimal  designs  are  shown  in  Fig.  14  .  The  topology 
designs  are  obtained  using  bilinear  4  node  quad,  B-spline  9  node  and  B-spline  16  node 
elements.  SIMP  interpolation  method  with  the  penalty  parameter  p  =  4  for  the  density 
function  and  the  allowable  material  volume  fraction  of  0.3  is  used.  It  can  be  observed 
that  with  the  use  of  sparse  mesh,  the  bilinear  quad  4-noded  elements  results  in  a  shape 
that  is  not  well  connected  and  have  problems  in  smooth  representation  of  the 
boundaries.  The  optimal  geometries  obtained  using  B-spline  elements  are  well 
connected  and  smoother  without  any  checkerboard  patterns. 
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(a)  (b)  (c) 


Fig.  15  Topology  results  for  a  mechanical  gripper  design  with  a  50  x  50  mesh  using  (a) 
Quad  4N  elements  (b)  B-spline  9N  elements  (c)  B-spline  16N  elements 

The  topology  results  with  an  increased  mesh  refinement  of  50  x  50  elements  are 
shown  in  Fig.  15  for  quadrilateral  4-node  (Q4),  B-spline  9-Node  and  B-spline  16-Node 
elements.  It  can  be  observed  that  with  the  increase  in  mesh  refinement,  even  the  Q4 
elements  converge  to  a  better  smooth  shape  and  the  B-spline  elements  also  converge  to 
better  smooth  shapes. 

To  evaluate  the  validity  of  the  designs  obtained  using  B-spline  elements  in  IBFEM, 
finite  element  models  similar  to  the  optimal  designs  were  created  using  the  commercial 
FEA  package  ABAQUS.  The  finite  element  model  of  the  optimal  design  along  with  the 
loads  and  boundary  conditions  are  shown  in  Fig.  16  .  Bilinear  quadrilateral  4-noded 
plane  stress  elements  are  used  for  the  analysis.  A  superimposed  image  of  the  deformed 
and  un-deformed  shapes  of  the  geometry  is  shown  in  Fig.  16  (b).  The  deformation  of  the 
mechanism  was  indeed  in  the  direction  as  intended.  Thus,  the  designs  obtained  for  the 
gripper  mechanism  are  indeed  valid. 
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Fig.  16  Results  from  ABAQUS  for  the  gripper  mechanism  (a)  FE  Model  of  the  gripper 
mechanism  with  loads  and  boundary  conditions  (b)  Deformed  shape  of  the  gripper 

mechanism 

Displacement  Inverter  Mechanism  Design 

The  feasible  domain  for  the  design  of  a  displacement  inverter  mechanism  is  shown 
in  Fig.  17  .  The  mechanism  is  supported  at  the  two  corners  along  its  left  edge  and  input 
forces  of  magnitude  50000  N  are  applied  at  the  middle  of  the  left  edge  as  shown  in  Fig. 
17  .  Displacements  are  expected  at  the  output  ports  at  the  middle  of  the  right  edge  in  the 
negative  x-direction.  Forces  of  magnitude  5000  N  are  applied  at  the  output  ports  in  the 
direction  opposite  to  the  direction  in  which  displacements  are  required. 


Fig.  17  Feasible  domain  for  displacement  inverter  design  with  a  30  x  30  mesh 

The  size  of  the  design  domain  is  5  x  5  m  as  shown  in  Fig.  17  .  Displacements  of 
magnitude  -1x10"^  (in  the  negative  x-direction)  are  specified  at  output  ports  and 
displacements  of  magnitude  1x10"^  are  specified  at  the  points  where  the  input  forces  are 
applied.  The  material  of  the  domain  is  assumed  to  be  steel  with  modulus  of  elasticity 
equal  to  200  GPa  and  the  Poisson's  ratio  of  0.3.  The  original  domain  has  been 
discretized  with  a  sparse  mesh  of  30  x  30  elements. 

The  topology  results  of  the  optimal  designs  for  the  inverter  mechanism  are  shown 
in  Fig.  18  .  The  topology  designs  are  obtained  using  bi-linear  4  node  quad,  B-spline  9 
node  and  B-spline  16  node  elements.  SIMP  interpolation  method  with  the  penalty 
parameter  p  =  4  for  the  density  function  and  the  allowable  material  volume  fraction  of 
0.2  is  used. 
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(a)  (b)  (c) 


Fig.  18  Topology  results  for  a  displacement  inverter  design  with  a  30  x  30  mesh  using  (a) 
Quad  4N  elements  (b)  B-spline  9N  elements  (c)  B-spline  16N  elements 


With  the  use  of  a  sparse  mesh  of  30  x  30  elements,  the  design  obtained  using 
bilinear  quad  4-noded  elements  is  not  well  connected  and  the  boundary  representation 
is  not  smooth.  On  the  other  hand,  B-spline  elements  result  in  geometries  that  have 
smooth  and  clear  boundaries.  Checkerboard  pattern  is  inherently  eliminated  in  B-spline 
elements. 

The  second  part  of  this  example  is  performed  on  a  similar  feasible  domain  with  a 
refined  mesh  discretization.  A  refined  mesh  discretization  of  50  x  50  elements  was  used 
to  validate  if  the  topologies  obtained  would  be  any  different  from  the  topologies 
obtained  using  a  sparse  mesh.  The  topology  results  of  the  optimal  designs  are  shown  in 
Fig.  19  .  With  a  refined  mesh,  the  design  obtained  using  bilinear  Quad  4-node  elements 
has  considerably  improved  with  clear  boundaries.  B-spline  elements  also  result  in 
geometries  that  have  smooth  and  clear  boundaries. 


(a)  (b)  (c) 


Fig.  19  Topology  results  for  a  displacement  inverter  design  with  a  50  x  50  mesh  using  (a) 
Quad  4N  elements  (b)  B-spline  9N  elements  (c)  B-spline  16N  elements 
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The  designs  obtained  using  B-spline  elements  in  IBFEM  were  validated  using  a 
commercial  FEA  package  to  evaluate  the  working  of  the  displacement  inverter 
mechanism.  Commercial  FEA  package  ABAQUS  is  used  to  evaluate  the  designs.  The 
finite  element  model  of  the  optimal  design  of  the  displacement  inverter  mechanism 
along  with  the  loads  and  boundary  conditions  is  shown  in  Fig.  20  (a).  Bilinear  4-noded 
quadrilateral  elements  are  used  for  the  analysis.  Fig.  20  (b)  shows  the  superimposed 
image  of  the  deformed  and  un-deformed  geometries.  The  tip  of  the  inverter  is  expected 
to  move  in  the  negative  x-direction  when  a  force  is  applied  in  the  positive  x-direction  on 
left  edge.  The  deformed  shape  shows  that  the  designs  obtained  for  the  inverter 
mechanism  are  valid. 

(a)  (b) 

Fig.  20  Results  from  ABAQUS  for  the  inverter  mechanism  (a)  FE  Model  of  the  inverter 
mechanism  with  loads  and  boundary  conditions  (b)  Deformed  shape  of  the  inverter 

mechanism 

B-spline  elements  thus  demonstrate  the  ability  to  obtain  the  optimal  shapes  even 
with  sparse  mesh  discretizations  when  compared  with  the  bilinear  quadrilateral 
elements  which  required  dense  mesh  discretization  to  obtain  similar  optimal  shapes. 

Flapping  Wing  Mechanism 

A  flapping  wing  mechanism  for  a  micro  air  vehicle  is  to  be  designed  to  obtain 
large  displacements  at  the  tip  of  the  wings.  The  wings  will  be  activated  by  an  magnetic 
actuator  placed  in  the  fuselage.  The  shape  of  the  casing  for  the  actuator  is  to  be  obtained 
using  topology  optimization  so  that  it  functions  as  a  compliant  mechanism  as  well  as 
the  support  for  the  wing.  The  feasible  domain  for  the  flapping  mechanism  with  a  mesh 
size  of  70  X  35  elements  is  shown  in  Fig.  21  .  A  polynomial  power  oip  =  3  and  a  volume 
fraction  of  0.5  is  used  to  obtain  the  optimum  topology  results.  The  material  of  the 
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structure  is  assumed  to  be  steel  with  a  modulus  of  elasticity  200  GPa  and  Poisson's  ratio 
of  0.3. 


Fig.  21  Feasible  domain  for  a  flapping  wing  mechanism  with  70  x  35  elements 

Displacements  of  magnitude  2x10"^  are  specified  at  the  wing  tips  and  the  entire 
structure  is  fixed  at  the  centre  of  the  bottom  edge.  The  topology  results  obtained  using 
bilinear  4-node  quadrilateral  elements,  B-spline  9N  elements  and  B-spline  16N  elements 
are  shown  in  Fig.  22  . 


Fig.  22  Topology  results  for  the  flapping  mechanism  for  a  70  x  35  size  mesh  and  a  volume 
fraction  of  0.5  (a)  Quad  4N  elements  (b)  B-spline  9N  elements  (c)  B-spline  16  N  elements 
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To  validate  the  design  of  the  mechanism  obtained  using  topology  optimization,  a 
finite  element  analysis  is  performed  on  the  final  topology  of  the  structure.  Fig.  23 
shows  a  superimposed  image  of  the  deformed  shape  on  the  optimum  structure.  As 
expected  a  displacement  of  2x10"^  is  obtained  at  the  wing  tips  as  shown  proving  the 
validity  of  the  design. 
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Fig.  23  Results  from  a  finite  element  analysis  on  the  optimum  structure  for  the  flapping 

mechanism 
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IV.  MAGNETOSTATIC  ANALYSIS  USING  IMPLICIT 
BOUNDARY  FINITE  ELEMENT  METHOD 

1.  Overview 

Magnetostatic  analysis  and  force  computation  for  magnetically  actuated  devices 
involves  modeling  an  assembly  of  components  with  different  material  properties.  When 
traditional  finite  element  method  is  used  for  such  analysis,  it  requires  a  conforming 
mesh  that  approximates  the  geometry  of  the  assembly.  The  mesh  must  contain  nodes 
along  the  external  boundaries  and  the  interfaces  between  parts.  The  edges  /  faces  of  the 
elements  must  approximate  these  boundaries  and  interfaces  as  shown  in  Fig.  24  .  Often 
the  geometry  is  not  well  approximated.  Generating  such  a  mesh  is  difficult  and,  despite 
decades  of  research,  3D  mesh  generation  (especially  using  hexahedral  elements)  is  still 
not  a  fully  automated  process  and  in  fact  requires  significant  user  input.  To  address 
mesh  generation  difficulties  several  meshless  methods  [30]  have  been  proposed  that  still 
need  a  well-placed  distribution  of  nodes  but  do  not  require  these  nodes  to  be  connected 
into  elements.  Some  of  these  methods  have  been  successfully  used  for  magnetostatic 
analysis  [31]-[36].  These  methods  use  interpolation  and  approximation  schemes  that  do 
not  need  connectivity  between  nodes.  However,  computationally  these  methods  are 
significantly  more  expensive  and  they  still  approximate  boundaries  and  interfaces  using 
nodes  along  them. 


Interface  boundaries 


(a)  33  elements  (b)  82  elements 


Fig.  24  2D  FEM  mesh 

An  alternate  approach  to  avoid  mesh  generation  difficulties  is  to  use  a  structured 
background  mesh  to  represent  the  solution  while  using  accurate  equations  of  curves 
and  surfaces  to  represent  the  boundaries.  A  structured  mesh  consists  of  uniform  regular 
shaped  elements  and  is  therefore  easy  to  generate.  Extended  finite  element  method  (X- 
FEM)  [37]- [39]  is  one  such  method,  which  uses  a  structured  mesh  and  implicit  equations 
for  the  boundaries  and  interfaces.  In  the  X-FEM  approach,  the  solution  is  enriched  near 
singularities  and  discontinuities  such  as  cracks.  An  important  application  of  this 
method  has  been  fracture  mechanics,  where  crack  propagation  [40]-[41]  is  simulated  by 
modifying  the  equations  of  the  crack  rather  than  regenerating  the  mesh.  Boundary  and 
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interface  conditions  have  been  imposed  using  Lagrange  multiplier  and  Penalty  methods 
for  X-FEM. 


The  Implicit  Boundary  FEM  (IBFEM),  uses  solution  structures  constructed  using 
implicit  equations  of  the  boundaries  to  enforce  boundary  and  interface  conditions.  This 
method  has  been  applied  to  2D  and  3D  elastostatics  and  steady  state  heat  transfer 
problems  [42] -[45].  Structured  mesh,  which  has  uniform,  undistorted  elements,  can  be 
used  for  the  analysis  because  the  implicit  boundary  method  does  not  require  nodes  on 
the  boundary  to  impose  boundary  conditions.  Structured  mesh,  such  as  the  examples 
shown  in  Fig.  25  is  easy  to  generate  since  all  elements  are  regular  shaped  and  the  grid 
does  not  have  to  conform  to  the  geometry. 


(a)  Inner  conductor :  Gridi 


(b)  Insulator :  Grid  2  (c)  Outer  conductor :  Grid  3  (d)  Coaxial  cable  :  Total  Grid 


(a)  2D  structured  mesh 


(A)  Material  1  (B)  Material  2  (C)  Multi-material 

(b)  3D  structured  mesh 

Fig.  25  Structure  mesh  for  multi-material  systems 


For  modeling  multiple  materials  and  assemblies,  a  separate  grid  is  generated  for 
each  material  or  part  as  shown  in  Fig.  25  .  Within  overlapping  elements  at  the  interface, 
the  piece-wise  interpolation  within  each  grid  is  combined  into  a  single  solution 
structure. 


2.  Governing  equation  and  Weak  form  for  2D  Magnetostatics 

Under  static  and  quasi-static  conditions,  the  governing  equation  for  2D 
magnetostatics  is 
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=  -J 


(4.1) 
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"1  da'" 
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j _ 

fi  ad^ 

^x^ 

dx^ 

where,  //  is  the  magnetic  permeability,  A  is  the  component  of  magnetic  vector 
potential  in  the  Xj  -direction  (the  direction  normal  to  the  plane  of  analysis)  and  J  is  the 
current  density  in  the  x, -direction.  In  the  finite  element  method,  essential  boundary 
conditions  are  specified  by  assigning  values  for  the  nodes  along  the  boundary.  When  a 
structured  grid  is  used  for  the  analysis,  there  may  not  be  any  nodes  available  on  the 
boundary.  Solution  structures  that  use  the  equation  of  the  boundary  to  impose  essential 
boundary  condition  have  been  used  by  several  authors  [46]-[48].  For  2D  magnetostatics, 
a  similar  solution  structure  for  the  X3  component  of  magnetic  vector  potential  could  be 
defined  as 

A(x)  =  D(x)A^  (x)  +  A‘‘  (x)  =  A^  (x)  +  A‘‘  (x)  (4.2) 

where,  4®  is  a  grid  variable  that  is  defined  by  piece- wise  interpolation  or  using  B- 
spline  approximation  [44]  over  a  structured  grid.  A‘'  is  the  boundary  value  function 
which  has  a  value  equal  to  the  prescribed  boundary  conditions  at  the  boundaries.  D(x) 
is  a  weighting  function  defined  such  that  D(x)  =  0  at  boundaries  where  essential 
boundary  conditions  are  applied  so  that  A=A‘'  at  these  boundaries.  The  boundary  value 
function,  A%  is  constructed  by  interpolating  nodal  values  within  elements.  The  nodal 
values  are  selected  such  that  at  the  boundary  it  has  a  value  equal  to  the  specified 
boundary  condition.  Note  that  Z)(x)  =  0  can  be  any  type  of  implicit  equation  of  the 
boundary  but  in  general  it  is  hard  to  construct  a  global  function  that  is  zero  only  at  the 
boundaries  with  essential  boundary  conditions.  Furthermore,  a  global  weighting 
function  can  lead  to  poor  convergence  especially  if  it  is  nonlinear.  In  the  implicit 
boundary  method,  approximate  step  functions  referred  to  as  Dirichlet  functions  or  D- 
functions  [43]-[44]  are  used  as  the  weighting  function.  At  any  given  point  xeRh  the  D- 
function  is  defined  as 


D(x)  = 


0 

l-(l-^(x)/Sf 

1 


^(x)  <  0 

0  <  <j)(x)  <  S 
<^(x)  >  5 


(4.3) 


where,  (f>(x)  =  Q  is  the  signed  distance  function  or  the  distance  from  the  boundary 
with  a  negative  value  if  the  point  x  is  outside  the  domain.  is  a  parameter  which 
controls  the  width  of  the  transition  band  over  which  the  D-function  transitions  from  0  to 
1.  In  the  limit  as  <^^0  the  D-function  approximates  the  Heaviside  step  function.  The 
advantage  of  using  an  approximate  step  function  as  weighting  function  is  that  it 
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transitions  from  0  to  1  within  boundary  elements.  Therefore,  they  can  be  locally  defined 
within  the  required  boundary  elements  and  for  all  other  elements  D(x)  =  1 .  This  implies 
that  the  internal  elements  are  not  influenced  by  the  weighting  function  and  are 
therefore  identical  to  traditional  finite  elements.  Moreover,  all  the  internal  elements  are 
identical  to  each  other  and  have  the  same  stiffness  matrix  since  they  have  the  same 
shape  and  size.  The  value  of  S  is  chosen  to  be  less  than  one-tenth  of  the  element  length 
in  our  numerical  implementation  so  that  the  D-function  closely  approximates  the 
Heaviside  step  function.  Using  the  solution  structure  defined  in  (2),  the  weak  form  of 
the  2D  magnetostatic  equation  can  be  derived  as 

^V(SA^)ju-W(A^)dV  = 

f  f  f  , 

j  (SA^)JdV  +  j  (SA^)H^SdS  -  J  V(SA^)ju~'V(A‘‘)dV 

V  V  V 


where,  SA'  is  the  virtual  magnetic  potential  vector  and  H,  is  the  tangential 
component  of  the  magnetic  field. 

The  grid  variable,  ,  is  interpolated  within  each  element  as  A^  ={n}^|a*|  where, 
{n}^  is  a  row  matrix  containing  the  shape  functions  and  |a*}  is  a  column  matrix 

containing  the  nodal  values  of  the  grid  variable.  Similarly,  the  boundary  value  function 
is  represented  within  each  element  as  yt"  ={n}^|a‘'|  where,  jA"}  is  column  matrix 

containing  the  nodal  values  assigned  such  that  A‘‘  has  the  prescribed  value  at  the 
boundary.  Note  that  the  same  shape  functions  are  used  to  interpolate  A^  and  A" .  If  all 
the  essential  boundary  conditions  are  homogeneous,  so  that  ^  =  0  is  the  only  prescribed 
boundary  condition,  then  the  boundary  value  function  A“  is  zero  everywhere  and  can 
be  eliminated  from  the  solution  structure.  Otherwise,  nodes  near  the  boundary  are 
assigned  values  of  A"  equal  to  the  prescribed  value.  For  2D  problems,  the  gradients  of 
the  boundary  value  function  A‘‘  is  expressed  as 


\dA^ 

[dx, 


dx^ 


'ax 


(4.5) 


The  gradients  of  the  homogenous  part  of  the  solution  A'  is  stated  as 


VA^  = 


dA^  M 
dx^  8x2 
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[B]  is  decomposed  into  two  matrices  [BJ  and  [BJ  such  that  only  [BJ  contains 
derivates  of  the  D-function  which  can  have  very  large  values  near  the  boundary. 


[BJ  =  .^.=D^ 

ox,  j 

^  and 


(4.7) 

(4.8) 


The  element  matrix  to  be  assembled  into  the  global  equations  can  be  defined  as 


[k^]  =  j  [b]'  [b]  =  [k[] + [k^  ] + [k;] 

[K[]=  j[B, ]"//-■  [B,]^Q, 

Q.e 

[k^]=  j[B, ]"//-■  [B,]^/Q, 

Q.e 

[k;  ]  =  j  ([B  J  //-■  [B,  ]  +  [B,  ]'  //-■  [B,  ]) 

Q-e 


(4.9) 

(4.10) 

(4.11) 

(4.12) 


[Bj]  which  contains  the  derivatives  of  D-function,  is  non-zero  only  within  the 

narrow  transition  band  near  the  boundary.  Therefore,  for  all  internal  elements  and 
boundary  elements  without  essential  boundary  conditions  [kj]  and  [kj]  are  zero.  For 

boundary  elements  [k[]  is  evaluated  by  subdividing  these  elements  into  triangles  and 

integrating  only  within  triangles  that  are  inside  the  geometry.  For  boundary  elements 
with  boundary  conditions,  the  volume  integral  for  computing  [k']  and  [Kj]  can  be 

converted  to  surface  integrals  because  they  contain  [b^]  which  is  non-zero  only  within 
the  narrow  transition  band  near  the  boundary.  The  components  of  [KJ]  can  be 
expressed  using  index  notation  as 


= 


where. 


"az)" 


|V^| 


dcj) 


(4.13) 


(4.14) 
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In  the  preceding  equation,  the  volume  integral  in  (11)  has  been  converted  into  a 
combination  of  surface  integral  along  the  boundary  and  an  integration  over  (p .  Similarly, 
the  components  of  [KJ]  can  be  stated  as 


^3,  =  JZ. 


dN, 


dx. 


dx. 


\dr^ 


(4.15) 


where. 


VO 


DD  1 


dcj) 


(4.16) 


All  components  of  [k;'],  [kj],  and  [Kj]  are  evaluated  using  Gaussian  quadrature. 

For  surface  integrals,  the  boundary  within  element  is  approximated  by  sufficiently 
small  straight  line  segments  to  achieve  accuracy. 


3.  Governing  equation  and  Weak  form  for  3D  magnetostatics 

Several  alternate  formulations  have  been  proposed  in  literature  for  3D 
magnetostatic  analysis  using  finite  element  method  [49]-[57].  A  formulation  based  on 
magnetic  vector  potential.  A,  was  used  in  our  implementation.  The  governing  equations 
for  3D  magnetostatics  expressed  in  terms  of  magnetic  vector  potential  is 

Vx(vVxA)  =  J  in  Q  (4.17) 

where,  Q  is  the  domain  of  analysis.  The  boundary  of  the  analysis  domain  T 
consists  of  regions  with  specified  natural  boundary  conditions  and  regions  that  are 
open  boundaries,  which  are  used  to  artificially  truncate  the  analysis  domain  when  in 
reality  it  extends  to  infinity.  Often  homogeneous  essential  boundary  conditions  are 
used  on  these  open  boundaries  as  an  approximation  if  the  boundary  is  far  away  from 
the  sources.  Several  special  techniques  for  modeling  such  open  boundaries  have  been 
developed  such  as  the  infinite  elements  and  asymptotic  boundary  condition  [58]. 
Natural  boundary  conditions  can  be  applied  on  boundaries  (denoted  as  T^  )  with 
known  tangential  component  of  the  magnetic  field  or  on  boundaries  (denoted  as  Tg ) 
with  known  normal  component  of  the  flux  density.  If  these  boundaries  are  planes  of 
symmetry  then  n  •  (V  x  A)  =  0  on  Tg  and  n  x  (vV  x  A)  =  0  T^ .  To  ensure  uniqueness  of 

the  solution,  the  following  essential  boundary  conditions  are  used  to  enforce  these 
conditions  [51]. 

nxA  =  0  onTg  (4.18) 
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n  •  A  =  0 


(4.19) 


onTjj 

Fig.  26  shows  an  example  domain  of  analysis  which  may  contain  regions  of 
different  materials  ( Qj  and  shown.  rj2  is  the  interface  surface  between  the  two 

sub-domains  Qj  and  Qj  shown  in  Fig.  26  .  At  the  interface,  the  tangential  component 
of  the  magnetic  field  and  the  normal  component  of  the  flux  density  are  continuous. 


Fig.  26  Analysis  domain  and  boundaries 

The  weak  form  for  these  governing  equations  and  boundary  conditions,  obtained 
using  the  weighted  residual  method  [59],  is 

|(Vx<JA)-(vVxA)JQ=  j  <JA(Hxii)/dr  +  |j^AJQ  (4.20) 

Q  r^+r^  Q. 

where,  SA  is  the  vector  weighting  functions.  This  weak  form  is  used  in  the 
traditional  FEM  to  compute  the  element  matrices  by  integrating  the  left  hand  side  over 
the  volume  of  each  element.  When  a  structured  mesh  is  used  for  the  analysis,  the 
boundaries  pass  through  the  elements  so  that  it  is  necessary  to  integrate  over  partial 
volume  of  the  element  that  is  inside  the  boundary.  Several  techniques  [60]-[61]  have 
been  developed  for  integrating  over  partial  elements  approximated  as  polygons. 
Alternatively,  the  partial  boundary  elements  can  be  subdivided  into  triangles  (for  2D) 
or  tetrahederons  (for  3D)  for  integration  purpose.  The  generated  triangles  or 
tetrahedrons  are  used  only  for  quadrature  and  not  to  represent  the  solution.  Even 
though  the  tessellation  of  the  boundary  elements  for  integration  approximates  the 
boundary,  the  size  of  the  triangles/  tetrahedrons  can  be  much  smaller  than  the  elements 
of  the  grid.  So  the  geometry  can  be  represented  reasonably  accurately  even  if  a  sparse 
mesh  is  used  for  the  analysis.  Essential  boundary  conditions  are  applied  in  traditional 
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FEM  by  assigning  values  to  the  nodes  on  the  boundary.  However,  when  a  structured 
mesh  is  used  there  may  not  be  nodes  available  on  the  boundary.  In  the  next  section,  a 
solution  structure  is  described  that  is  constructed  using  the  implicit  equation  of  the 
boundary  to  enforce  essential  boundary  conditions. 

Several  structured  mesh  based  approaches  for  analysis  have  used  solution 
structures  [42]-[48]  constructed  using  implicit  equations  of  boundary,  to  impose 
essential  boundary  conditions.  Here  we  present  the  implicit  boundary  method,  where 
step  functions  are  used  as  implicit  equations  to  construct  solution  structures.  For  three- 
dimensional  magnetostatics,  a  solution  structure  for  the  magnetic  vector  potential  A(x) 
could  be  defined  as 

A(x)  =  [D(x)]  A^(x)  +  A‘'(x)  =  A^(x)  +  A“(x)  (4.21) 

In  the  preceding  equation,  A^  is  a  grid  variable  vector  that  is  defined  by  piece- 
wise  interpolation  or  using  B-spline  approximation  [44]  over  a  structured  mesh  and 
D(x)  is  diagonal  matrix  whose  components  are  defined  such  that  D;;(x)  =  0  at  the 
boundaries  where  essential  boundary  conditions  are  applied  on  the  ith  component  of 
A  .  This  ensures  that  A.  =  at  the  boundary.  A*  is  the  homogenous  part  of  the 

solution  and  A''  is  the  boundary  value  function.  The  boundary  value  function  is 
defined  such  that  it  has  value  equal  to  the  prescribed  boundary  conditions  at  the 
boundaries.  Therefore,  diagonal  components  of  the  D  matrix,  D|.(x) ,  are  implicit 
equations  of  the  boundaries  on  which  essential  boundary  conditions  are  applied.  It  is 
hard  to  construct  such  a  function  that  is  zero  only  at  the  boundaries  with  essential 
boundary  conditions.  Moreover,  if  this  weighting  function  is  nonlinear  and  defined 
globally,  then  this  solution  structure  can  lead  to  poor  convergence.  In  the  implicit 
boundary  finite  element  method,  approximate  step  functions  are  used  as  the  implicit 
equation.  At  any  given  point  x  g  R’  this  step  function,  referred  to  as  the  D-function,  is 
defined  as 


0 

l-(l-^(x)/^f 

1 


(^(x)<0 
0  <  <^(x)  <  5 
(l){x)  >  5 


(4.22) 


where,  is  the  signed  distance  function  for  the  boundary.  The  signed  distance 
function  for  a  boundary  is  evaluated  at  any  point  x  as  the  distance  of  the  point  from  the 
boundary.  The  function  has  a  negative  value  if  the  point  is  outside  the  domain.  The  step 
function  D;;(x)  transitions  from  0  to  1  over  a  band  whose  width  is  controlled  by  the 
parameter  S .  In  the  limit  as  ^  0  the  D-function  approximates  the  Heaviside  step 
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function.  In  our  numerical  implementation,  the  value  ot  the  parameter  5  is  chosen  to  be 
less  than  one-tenth  of  the  element  length.  The  advantage  of  using  an  approximate  step 
function  as  weighting  function  is  that  it  transitions  from  0  to  1  within  boundary 
elements.  Therefore,  they  can  be  locally  defined  within  these  elements.  For  all  other 
elements,  which  do  not  have  a  boundary  with  prescribed  essential  boundary  conditions, 
we  can  set  D;;(x)  =  1.  This  implies  that  the  internal  elements  are  not  influenced  by  the 
weighting  function  and  since  the  structured  mesh  is  made  of  uniform  elements  that  are 
identical  to  each  other,  all  the  internal  elements  have  the  same  stiffness  matrix. 
Substituting  the  solution  structure  (5)  into  the  weak  form  (4),  a  modified  weak  form  of 
the  3D  magnetostatic  equation  can  be  derived  as 

I  (  V  X  dA* )  •  (vV  X  A*  )dQ  =  I  JSA^dQ  - 1  (  V  x  dA* )  •  (vV  x  A“  (4.23) 

Q.  no. 


where,  dA*  is  the  virtual  magnetic  potential  vector. 

The  grid  variable  vector,  A^  ,  is  interpolated  within  each  element  as 
a''=K{  A^l  where,  [Nf  is  a  matrix  containing  the  shape  functions  and  is  a 
column  matrix  containing  the  nodal  values  of  the  grid  variable  vector.  For  brick 
elements  with  8  nodes,  the  size  ot  {a^|  is  24  because  the  nodal  degree  of  freedom  is  3. 

Similarly,  the  boundary  value  function.  A",  is  defined  by  interpolating  nodal  values 
within  each  element  as  A”=[Nf{  A"  I  where,  {A"!  is  column  matrix  containing  the 

nodal  values  of  A" .  These  nodal  values  are  assigned  such  that,  at  the  boundaries,  this 
function  will  have  a  values  prescribed  by  the  boundary  condition.  Using  the  solution 
structure,  the  magnetic  flux  for  the  boundary  value  function  can  be  derived  as 

Vx  A"  =[B'^]|A‘'}  (4.24) 


The  'curl'  matrix  [b'^]  for  the  boundary  value  function 


[B 

/]  ... 

[«/]]. 

where. 
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dz 
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dN, 
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dN, 
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dN, 
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0 

for  ieZ  =  [l,  A] .  The  curl  ot  A*  can  be  computed  as 


is  defined  as 


(4.25) 


45 


VxA*  =[B^]{a^} 


(4.26) 


For  convenience,  [b*^]  is  defined  as  a  sum  of  two  matrices  such  that  the  first  one 

only  contains  derivatives  of  the  shape  function  and  the  second  matrix  contains  the 
derivatives  of  the  D-function.  [b^^]  =  [b,‘^]+[b/]  ,  where 


8.q=[ 

'[B^q  [ 

B/]  ...  1 

;B,.q] 

v]-[ 
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0 

(4.27) 


(4.28) 


(4.29) 


In  the  preceding  equations,  i  =  ,  where,  n  is  the  number  of  nodes  per 

element.  The  element  matrix  that  is  assembled  into  the  global  equations  can  be  defined 
as 


[k-  ]  =  J  {[b^V  [sq)  rfn,  =  [k;  ] + [k-  ] + [k-  ] 

[Kq=f{[B/7K[B/]jrfn, 

[K;]=j({[Bqv[B,q)q[B/j.[Bq))rfn, 


(4.30) 

(4.31) 

(4.32) 

(4.33) 
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Since  [b/]  contains  the  derivates  of  D..(x),  it  is  non-zero  only  within  the  narrow 

transition  band  near  the  boundary.  Therefore,  for  all  internal  elements  and  boundary 
elements  without  essential  boundary  conditions  [Kj]  and  [k']  are  zero.  Within  the 

transition  band,  the  derivatives  of  D(x)  can  have  large  magnitude.  For  the  boundary 
elements  with  boundary  conditions,  the  volume  integral  for  computing  [k']  and  [k*] 
must  be  converted  to  surface  integrals  as  follows  to  compute  them  accurately. 

I  ''[»=" (4.34) 

=  (4.35) 

To  derive  the  preceding  equations,  we  make  use  of  the  fact  that  is  zero 
except  in  the  narrow  band  0  <  J .  Therefore,  the  volume  integral  is  converted  into  a 
surface  integral  along  the  boundary  and  an  integral  over  the  transition  band  (normal 
to  the  surface).  Note  that  if  (z>  is  a  signed  distance  function  then  |v^|  =  l.  If  the  width  of 
the  band  S  is  very  small,  then  one  can  assume  that  the  shape  functions  are  constant 
within  the  band,  allowing  the  integral  over  (j)  to  be  determined  analytically. 
Alternatively,  the  integration  over  (zi  can  also  be  evaluated  numerically. 

4.  Solution  structure  for  multi-material  models 

At  the  interface  between  materials  with  different  magnetic  permeability,  the 
normal  component  of  magnetic  field  and  the  tangential  component  of  flux  density  can 
be  discontinuous.  The  tangential  component  of  the  magnetic  field  and  normal 
component  of  the  flux  density  are  continuous.  The  required  interface  conditions, 
expressed  in  terms  of  the  magnetic  vector  potential  are 

nx(vjVx  Aj)  =  Aj)  (4.36) 

n-(Vx  Aj)  =  n-(Vx  A2)  (4.37) 

It  is  obvious  that  if  the  magnetic  vector  potential  is  continuous,  that  is  Aj  =  A2 , 
then  the  second  condition  (21)  is  automatically  satisfied.  The  requirement  for  continuity 
of  the  tangential  component  of  the  magnetic  field  (20)  requires  that  the  derivatives  of 
the  vector  potential  should  be  discontinuous.  To  allow  this  discontinuity  in  the 
magnetic  field,  separate  grids  are  used  for  each  material  as  shown  in  Fig.  25  .  At  the 
interface,  the  elements  from  neighboring  grids  overlap.  A  solution  structure  is  needed 
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for  overlapping  elements,  which  ensures  that  the  vector  potential.  A,  is  continuous 
while  flux  density  B  and  magnetic  field  H  can  be  discontinuous.  The  following  solution 
structure  was  used  for  interface  elements 

=  (1  -  D(x))  +  D(x)  A"""  (4.38) 


where.  A®'  is  the  field  interpolated  or  approximated  within  the  element  from  grid 
i,  (i  =  1,  2),  D(^^j(x))  is  the  approximate  step  function  defined  in  (6)  and  <f)^{x)  is  the 

implicit  equation  of  the  interface  curve  (represented  using  signed  distance  function). 
Similar  solution  structure  has  been  used  to  model  material  discontinuity  for  elasticity 
problems  [45].  The  solution  structure  in  (22)  blends  the  solutions  from  the  two  grids 
such  that  the  vector  potential  is  continuous  at  the  interface.  Note  that  this  method  for 
blending  the  solutions  uses  a  partition  of  unity  as  weighting  functions.  This  solution 
structure  ensures  the  continuity  of  the  solution  throughout  the  analysis  domain.  It  also 
allows  the  derivatives  (and  magnetic  field  and  flux  density)  to  be  discontinuous  at  the 
interface.  The  gradient  of  the  vector  potential  components  are 
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(4.39) 


For  small  values  of  5,  the  gradient  of  D(x)  can  have  large  magnitude  which  in 
turn  acts  as  a  penalty  that  tries  to  enforces  A*‘  =  at  the  interface  so  that  the  second 
and  third  terms  cancel  each  other.  Therefore,  the  grid  variable  is  approximately 
continuous  while  the  vector  potential  is  exactly  continuous.  The  slope  of  the  grid 
variable  and  therefore  the  vector  potential  can  be  continuous  or  discontinuous  as 
dictated  by  the  equilibrium  equations. 

Substituting  (22)  into  the  weak  form,  the  element  matrices  for  the  interface 
elements  can  be  computed.  Again,  since  gradient  of  D(x)  is  very  large  near  the 
interface,  it  helps  to  decompose  the  matrices  into  terms  that  only  contain  derivatives  of 
the  shape  functions  and  those  that  contain  derivative  of  D(x) .  As  illustrated  in  the 
previous  section,  all  the  terms  that  involve  derivates  of  D(x)  can  be  converted  from 
volume  integrals  into  surface  integrals  for  accurate  numerical  evaluation.  These 
techniques  are  described  in  detail  in  [45]  tor  elastostatics  and  have  been  adopted  here 
for  magnetostatics. 


5.  Magnetic  force  computation 

Several  techniques  for  computing  magnetic  forces  can  be  found  in  literature  [62]. 
These  include  the  Maxwell's  stress  tensor,  equivalent  source  method  and  the  virtual 
work  principle  [63]-[66]  to  list  a  tew.  These  approaches  have  been  implemented  using 
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FEM  and  therefore  can  also  be  used  with  the  implicit  boundary  FEM.  Since  we  assume 
that  the  exact  equation  of  the  surface  is  available  (preferably  as  a  parametric  equation), 
it  is  easier  to  implement  a  method  that  integrates  surface  force  densities  to  compute  the 
nodal  forces.  The  basic  equation  for  magnetic  force  density  [67]  can  be  derived  from 
energy  balance  equations  and  is  given  by 

K=-\h-vh+}>,'b 

^  (4.40) 

The  first  term  in  (19)  involves  the  gradient  of  magnetic  permeability.  Therefore, 
this  term  is  of  significance  only  at  the  boundary  between  ferromagnetic  materials  and 
surrounding  non-ferromagnetic  materials.  The  second  term  is  a  body  force  that  exists 
on  current  carrying  conductors  in  the  presence  of  magnetic  flux  density.  To  compute 
the  structural  response  due  to  these  magnetic  forces,  a  subsequent  solid  mechanics 
analysis  is  necessary.  The  weak  form  for  solid  mechanics  problems  is  the  principle  of 
virtual  work  which  can  be  stated  as  follows: 

{8£}^  [C]  {£}  dV  =  t  •  6ndS  (4.41) 

The  first  term  on  the  right  hand  side  of  the  weak  form  is  the  virtual  work  done  by 
magnetic  forces.  This  term  can  be  evaluated  as  follows: 

4  •  8uJF  =  jubwdV  + ( J  X B) •  8udF  (4.42) 

If  we  assume  that  the  permeability  is  constant  within  the  materials  then  the  first 
term  in  (21)  makes  a  contribution  only  at  the  boundary.  For  a  ferromagnetic  object  with 
permeability,  |a, ,  surrounded  by  a  medium  whose  permeability  is  ,  the  permeability 
can  be  considered  to  change  from  one  value  to  the  other  over  a  band  along  the 
boundary  whose  width,  measured  in  the  normal  direction,  is  A« .  The  limit  of  A«  0 
represents  the  discontinuous  variation  at  the  boundary.  The  gradient  of  the 
permeability  within  this  band  can  then  be  written  as: 

rr  - 

Vp  =  — n  (4.43) 

dn 

The  unit  vector  n  is  the  direction  normal  to  the  boundary  between  the  two 
materials  with  different  permeability  and  points  in  direction  of  increasing  permeability. 
Using  this  expression  in  the  first  term  of  (21),  we  get. 
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Therefore  the  volume  integral  can  be  converted  to  a  surface  integral  where,  is 
the  surface  separating  two  materials  with  different  permeability.  The  virtual  work  due 
to  magnetic  body  forces  can  now  be  evaluated  as 

I  •  8uJF  =  •  bwdV  +  f,  •  bndV  (4.45) 

where,  T  =  JxB  is  body  force  and  =  is  a  surface  force  density  (or 

traction).  The  surface  force  density  term  can  be  evaluated  by  expressing  the  square  of 
magnetic  field  as  a  function  of  permeability.  If  there  is  no  surface  current  then  the 
tangential  component  of  the  magnetic  field  does  not  vary  across  the  boundary  and  can 
be  treated  as  a  constant.  Similarly,  the  normal  component  of  the  magnetic  flux  density  is 
constant  (by  Gauss's  law)  and  does  not  vary  across  the  boundary  even  though  the 
permeability  is  different  on  the  two  sides  of  the  boundary.  The  square  of  the  magnetic 
field  can  be  expressed  as  the  sum  of  the  squares  of  the  tangential  and  normal 
components  of  the  field  as: 

=  Hf  +  Hi  =Hl+^  (4.46) 

h 


Both  the  tangential  component  of  magnetic  field  ( H,  )  and  the  normal  component 
of  the  magnetic  flux  density  ( )  can  be  treated  as  constants  for  the  integration  in 
computing  surface  traction  since  these  quantities  do  not  vary  in  the  direction  normal  to 
the  interface  between  the  two  materials.  Using  the  preceding  equation  for  the  square  of 
magnetic  field  to  compute  the  surface  traction  due  to  magnetic  forces  we  get 
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where,  and  n,  are  the  permeability  and  the  outward  normal  vector  of  the  ith 
material  and  is  the  tangential  component  of  the  magnetic  flux  density  within  the  ith 
material.  If  it  follows  that  /,  <  0  therefore  the  direction  of  the  surface  traction  is 

opposite  to  n,  which  means  that  it  acts  in  the  direction  of  decreasing  permeability.  In 
[66]  an  expression  for  the  surface  force  densities,  similar  to  (26),  between  two  linear 
media  has  been  deduced  from  a  more  general  expression  for  magnetic  force. 

At  the  interface  elements,  magnetic  forces  can  be  computed  by  integrating  the 
magnetic  force  density  over  the  interface  boundary.  The  nodal  forces  at  the  boundary 
elements  of  each  grid  can  be  computed  by  integrating  over  the  piece  of  the  boundary 
that  passes  through  the  element.  In  other  words,  for  a  boundary  element  whose 
material  property  is  //,  the  nodal  forces  are  computed  as: 
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(4.48) 


The  unit  normal  n,  +  is  constructed  such  that  it  points  outwards  from  the 
material  or  part  boundary.  The  tangential  and  normal  components  of  the  flux  density 
are  computed  from  the  vector  potential  as 


dA  dA 

Bi,=—na+—na 
0X2  ox^ 

dA  dA 

= - n,. - n,^ 

ox.  ox, 


(4.49) 


In  traditional  FEM  the  integration  over  the  interface  boundary  requires  finding  the 
elements  at  the  interface  and  determining  which  edge  or  face  lies  along  the  interface.  In 
the  implicit  boundary  approach,  since  the  equations  of  the  interface  are  available  for 
each  part/material,  it  is  very  easy  to  integrate  (27)  over  these  boundaries  to  evaluate  the 
nodal  forces  for  each  part  separately.  The  boundary  passing  through  each  element  is 
approximated  by  straight  lines  or  triangles  for  the  purpose  of  integration  and  Gauss 
Quadrature  is  used  to  perform  the  integration. 


6.  Results  and  discussion 

Several  examples  are  presented  here  that  were  used  to  validate  the  implicit 
boundary  method.  Some  of  these  examples  have  analytical  solutions  with  which  to 
compare  the  computed  results.  This  allows  us  to  not  only  verify  the  accuracy  of  the 
results  but  also  study  the  rate  of  convergence  and  compare  with  similar  results  using 
the  traditional  finite  element  method.  Some  examples  are  modeled  using  2D  and  3D 
elements  for  comparison  and  validation. 
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Coaxial  cable 


A  coaxial  cable,  which  consists  of  an  inner  conductor,  an  insulator,  and  an  outer 
conductor,  is  modeled  as  shown  in  Fig.  27  using  three  separate  structured  grids.  Due  to 
circular  symmetry  of  the  geometry,  only  a  quadrant  of  the  coaxial  cable  cross-section  is 
created.  The  radii  of  the  inner  conductor,  the  insulator  and  the  outer  conductor  are  a ,  b 
and  c .  The  inner  and  outer  conductors  carry  the  same  amount  of  total  current  in 
opposite  directions.  The  total  current  flowing  through  each  conductor  is  /  .  The  current 
flows  in  the  axial  direction  (z-direction)  and  the  current  density  is  assumed  to  be 
uniform. 


N 

\ 

\ 

\ 

V 

■ 

\ 

T 

1 

L 

1 

(a)  Inner  conductor :  Grid1  (b)  Insulator :  Grid  2 

(c)  Outer  conductor : 

:  Grid  3  (d)  Coaxial  cable  :  Total  Grid 

Fig.  27 

structured  mesh 

The  analytical  solution  of  the  magnetic  field  in 

circumferential  direction  can  be 

derived  as 

r(27ra^^  I 

,()<r  <a 

[iTrr)  *  I 

,a<r  <b 

{c^  {2nr)^  1 

,b<r  <c 

(4.50) 
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VI 

where,  r  is  the  radial  distance.  The  following  values  of  current  and  radii  were 
used  in  the  numerical  model:  /  =  lOOOd ,  a  =  0.5 ,  b  =  \,  and  c  =  1 .5  mm. 

Fig.  28  shows  the  magnitude  of  the  magnetic  field  that  was  computed  using  the 
quadratic  B-spline  elements.  It  shows  that  the  maximum  magnetic  field  value  is  at  the 
interface  between  the  inner  conductor  and  the  insulator  and  has  a  value  of  3.182x10^ 
A/mm.  This  is  very  close  to  the  value  obtained  from  the  analytical  solution  which  is 
3.183x10^  A/mm. 
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Fig.  28  Magnitude  of  the  magnetic  fieid 

For  3D  analysis  it  is  necessary  to  first  compute  the  current  density  by  solving  the 
electrostatics  equations.  Then  using  the  computed  current  density,  magnetic  field  is 
obtained  through  3D  magnetostatic  analysis.  The  governing  equations  for  the  electro¬ 
magnetostatic  problem  is 

-V-(aVF)  =  0  in  Q 

(31) 

Vx(vVxA)  =  J  in  Q 


where,  the  current  density  for  magnetostatics  is  J  =  -oVV .  Fig.  29  shows  the 
coaxial  cable  model  using  three  separate  structured  grids. 


A.  Inner  conductor  B.  Insulator 


C.  Outer  conductor 


Fig.  29  3D  coaxial  cable  model  with  the  structured  grid 

The  electric  conductivity  in  the  conductors  is  set  to  \0^S/mm  and  lO~^S/mm  in 
the  insulator.  In  order  to  obtain  7  =  1000^,  the  voltage  difference  in  the  top  and  the 
bottom  surfaces  is  set  to  0.25  V  in  the  inner  conductor  and  0.05  in  the  outer  conductor. 
The  current  density  of  the  inner  conductor  is  computed  as  1273  A/mm2  and  254.6 
A/mm2  in  the  outer  conductor  to  carry  the  amount  of  current. 
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Fig.  30  Current  density  in  z-direction  for  3D  coaxial  cable 


Fig.  30  shows  computed  current  density  in  z  direction  by  3D  electrostatic  analysis. 
The  calculated  current  densities  of  the  inner  conductor  and  the  outer  conductor  are  1259 
and  365  A/ 
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Fig.  31  Magnetic  field  for  3D  coaxiai  cabie 

Fig.  31  shows  the  magnitude  of  the  magnetic  field  that  was  computed  using  8 
brick  node  elements.  It  shows  that  the  maximum  magnetic  field  value  is  at  the  interface 
between  the  inner  conductor  and  the  insulator  and  has  a  value  of  3.35lxio^^/mm  .  This  is 
close  to  the  value  obtained  from  the  analytical  solution  which  is  3.183xio^.4/wm . 
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Fig.  32  Magnetic  field  versus  radius 


Fig.  32  shows  the  magnetic  field  in  the  hoop  direction  varying  with  the  radius.  In 
the  figure,  results  obtained  by  hexahedral  8-node  elements  are  denoted  as  H8.  This 
element  provides  a  piece-wise  trilinear  interpolation  for  the  vector  potential.  Therefore, 
the  derivatives  are  not  continuous  without  smoothing.  The  same  results  after 
smoothing  are  denoted  in  the  figure  as  H8S.  After  smoothing,  the  result  of  IBFEM  is 
very  close  to  the  analytical  solution.  Fig.  33  shows  the  convergence  of  HI  error  norm  for 
this  problem  using  trilinear  (H8)  elements  which  is  the  root  mean  square  error  in  flux 
density  field  over  the  domain  and  can  be  defined  as 


(33) 


where,  B^^  is  the  exact  value  of  the  magnetic  flux  from  the  analytical  solution  and 
B*  is  the  corresponding  computed  value. 
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Fig.  33  Convergence  of  H1  norm6.3.  Plunger  solenoid  actuator 
Switched  Reluctance  Motor 

A  2D  planar  model  of  the  Switched  Reluctance  Motor  (SRM)  is  shown  in  Fig.  34  . 
SRM  is  a  DC  motor  where  the  stator  has  windings  around  the  poles  while  the  rotor  does 
not  have  any  windings.  Current  is  applied  to  the  coils  around  the  poles  of  the  stator 
sequentially  to  produce  a  torque  on  the  rotor  as  it  rotates. 


(a)  Aligned  Position  (b)  Unaligned  Position 
Fig.  34  Switched  reluctance  motor 

Reference  [68]  provides  the  dimensions  of  the  motor  that  are  used  in  this  example. 
The  stator  and  the  rotor  are  made  of  iron  with  relative  permeability  of  2000.  A  quarter 
of  the  motor  is  modeled  in  both  its  aligned  and  unaligned  orientation.  The  number  of 
turns  in  the  coil  is  assumed  to  be  500  and  the  current  is  2  A.  Currents  only  flows  into 
coils  attached  to  the  top  pole  of  the  stator.  Essential  boundary  condition  (A=0)  is 
imposed  on  peripheries  of  the  shaft  and  the  stator  and  the  vertical  axis.  B-spline 
elements  are  used  for  the  analysis  so  that  accurate  results  are  obtained  even  with  a 
sparse  grid  as  shown  in  Fig.  35  .  The  flux  lines  computed  here  are  similar  those 
computed  by  FEM  [68].  Note  that  the  air  gap  is  very  small  compared  to  the  average 
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element  size  in  the  grid.  Despite  the  geometric  complexity  and  large  number  of  parts 
and  materials  involved,  this  approach  for  analysis  yields  results  using  structured  grids 
comparable  to  results  from  traditional  FEM  using  conforming  mesh. 

v<><(9r 

v^tur  rai«i  I  ak-  Ptttr 


IJ.7STE-2 

J.4UE-2 

J.03E  2 

I2.2W£-a 

1.8«E-2 

|EJ»7E-2 

EllJE  2 

}"■ 


il5E-S 

iA99tLi 

E282E-3 

LMSEJ 

i{t.483e-4 

8.314E-4 

4.t44E-4 

I1,??5E4 

-1.3I46E-5 


Fig.  35  Field  lines 


Iron  block  in  a  homogeneous  magnetic  field 

The  example  of  an  iron  cube  in  air  subject  to  homogenous  magnetic  field  has  been 
used  in  literature  to  verify  a  variety  of  formulations.  Fig.  36  shows  one-eighth  of  the 
system  modeled  considering  its  symmetry.  The  relative  permeability  of  iron  cube  is 
1000.  The  modeled  region  is  subjected  to  a  homogenous  magnetic  flux  density  in  the 
z-direction. 
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Fig.  36 


Iron  cube  in  homogeneous  magnetic  field 


The  half-length  of  the  iron  cube  edge  is  'a'.  The  symmetry  plane  are  x=0,  y=0  and 
z=0.  The  planes  x=b,  y=b,  and  z=b  represent  the  far  boundaries.  A  homogeneous 
magnetic  field  is  applied  in  the  z-direction  with  the  aid  of  boundary  conditions.  On  the 

B  b 

far  boundaries,  the  Dirichlet  boundary  conditions  are:  and  A^  =0  on  x=b  and 

B  b 

A  =  — ^  and  A  =  0  on  v=b.  The  dimensions  used  are  a  =  20mm ,  b  =  40mm  and  the 

X  2  " 

flux  density  magnitude  is  5o  =  1.0T  .  In  addition  to  this,  the  following  essential 
boundary  conditions  are  also  imposed 


nx  A  =  0 

on  Tg  (x=0  and  y=0) 

(29) 

o 

II 

on  Tjj  (z=0  and  z=b) 

(30) 

As  the  magnetic  flux  density  only  exists  in  the  z-direction,  the  normal  components 
of  magnetic  flux  density  must  be  zero  on  the  symmetry  planes  and  the  tangential 
component  of  magnetic  field  must  be  zero  on  planes  normal  to  the  z  direction. 
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A.Cube 


B.  Octagonal  Prism 


C.  Cylinder 


Fig.  37  Iron  objects  with  the  same  grid  density 


In  addition  to  modeling  the  iron  cube,  we  have  modeled  two  other  shapes  as 
shown  in  Fig.  37 ,  where  the  iron  part  is  modeled  as  an  octagonal  prism  and  a  cylinder. 
The  height  of  the  parts  is  the  same  and  equal  to  20mm.  The  edge  length  of  the  right 
octagonal  prism  is  20mm  and  the  cylinder  has  radius  =  20mm.  The  total  number  of 
elements  used  in  the  model  is  12167. 


A.  Cube 


B.  Octagonal  Prism  C.  Cylinder 


Fig.  38  Cross-sections  with  the  line  y=z=10mm 


As  shown  in  Fig.  38  ,  along  the  line  y=z=10mm,  the  interface  between  iron  and  air 
is  at  x=20mm  for  the  cube,  at  x=24.14mm  for  the  right  octagonal  prism  and  at 
x=17.32mm  for  the  cylinder.  Using  the  same  number  of  elements  and  the  same 
boundary  conditions  for  all  three  cases,  the  variation  on  the  three  components  of 
magnetic  flux  B  along  the  line  y=z=10mm  are  obtained  as  shown  in  Fig.  39 
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Fig.  39  Components  of  B  along  the  line  y=z=10mm.  (a)  Bx  (b)  By  (c)  Bz 

Fig.  39  (a)  shows  that  only  the  magnetic  flux  density  in  x-direction  is  continuous 
when  the  shape  of  the  iron  is  cubic  because  the  normal  direction  of  the  interface  is  along 
the  x-direction.  The  other  components  are  tangential  components  and  are  discontinuous 
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as  expected.  Similarly,  for  other  shapes  none  of  the  components  of  B  are  normal 
components  and  therefore  none  of  them  are  continuous  at  the  interface.  The  results 
obtained  for  this  example  are  similar  to  that  obtained  using  the  traditional  FEM  that 
uses  conforming  mesh. 

Plunger  solenoid  actuator 

Solenoid  actuators  are  designed  to  produce  small  linear  motion  of  an  armature 
and  are  of  several  types  depending  on  the  shape  of  the  armature.  A  solenoid  actuator 
with  a  plunger  armature  [69],  as  shown  in  Fig.  40  is  considered  in  this  example.  The 
armature  and  stator  are  made  of  steel  laminates  in  order  to  reduce  eddy  current  effect. 
The  stator  has  solenoidal  coil,  wound  in  the  shape  of  a  cylinder,  or  parallelepiped. 


Fig.  40  Plunger  solenoid  actuator 

In  Fig.  41  (a)  the  plunger  is  cylindrical  and  the  solenoid  is  axisymmetric,  while  in 
Fig.  41  (b),  the  plunger  has  square  cross-section  with  rounded  edges.  We  have  referred 
to  the  later  actuator  with  square  cross-section  as  the  brick  plunger  actuator.  For  both 
models,  the  number  of  turns  N=400  and  the  current  I=4A.  The  relative  permeability  of 
the  stator,  the  armature  and  the  stopper  is  =  2000  and  //^  =  1  in  the  coil. 
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A.  Cylindrical  plunger 


B.  Bricked  plunger 


Fig.  41  Top  view  of  piunger  actuators 


A.  Cylindrical  Plunger  B.  Bricked  Plunger 

Fig.  42  3D  soiid  modeis  of  the  soienoid  actuators 

Fig.  42  shows  3D  solid  models  of  the  solenoid  actuator  model,  which  were  used 
for  the  analysis.  In  order  to  reduce  the  computational  effort,  only  a  fourth  of  the  whole 
system  is  modeled.  The  current  density  in  the  coil  is  in  the  circumferential  direction  and 
is  obtained  by  first  solving  an  electrostatic  problem.  A  voltage  difference  of  0.0546  V 
was  applied  between  the  symmetry  planes  of  the  coil.  Assuming  the  conductivity  of  the 
coil  to  be  <j  =  \(f  S ! m  the  current  density  obtained  will  be  on  average  the  same  as  the 
current  density  applied  for  the  2D  axisymmetric  model.  The  same  voltage  boundary 
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conditions  are  applied  for  the  brick  plunger  actuator.  Using  3D  electrostatic  analysis, 
the  computed  current  density  is  obtained  as  shown  in  Fig.  43  . 
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A.  Cylindrical  Plunger 


B.  Bricked  Plunger 


Fig.  43  Magnitude  of  current  density 


The  computed  current  density  is  then  used  to  perform  the  magnetostatic  analysis, 
to  compute  magnetic  flux  density  and  magnetic  field.  At  the  symmetric  planes,  the 
essential  boundary  conditions,  nxA  =  0  ,  are  applied.  The  computed  magnetic  flux 
density  is  shown  in  Fig.  44  . 
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A.  Cylindrical  Plunger 


B.  Bricked  Plunger 


Fig.  44  Magnetic  flux  density  in  y  direction 


For  the  axisymmetric  case,  using  the  reluctance  method,  the  magnetic  flux  density 
in  the  air  gap  can  be  calculated  as  B=0.1715  T  and  the  magnetic  force  is  F  =14.7  N.  Brauer 
[69]  has  also  provided  2D  FEM  results;  B=0.170  T  and  F  =19.34  N.  The  magnetic  flux 
density  computed  using  IBFEM  is  also  roughly  in  this  range.  The  computed  magnetic 
field  density  in  the  y-  direction  is  shown  in  Fig.  45  . 
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A.  Cylindrical  Plunger 


B.  Bricked  Plunger 


Fig.  45  Magnetic  field  in  y  direction 


Fig.  45  (A)  shows  that  the  computed  magnetic  field  in  y  direction  is  close  to  the 
value  of  1.364x10^  A/m  computed  using  the  reluctance  method  [69].  Using  the  magnetic 
field  and  flux  density,  the  computed  magnetic  force  on  the  cylindrical  armature  is  18.33 
N,  which  is  quite  close  to  the  force  calculated  by  2D  FEM.  For  the  brick  plunger 
armature,  the  computed  force  is  19.56  N.  According  to  this  analysis,  the  force  on  the 
brick  shaped  armature  is  larger  than  on  the  cylindrical  armature  because  the  cross- 
sectional  area  of  the  bricked  armature  (  1.592x10"^  m^  )  is  larger  than  that  of  the 
cylindrical  armature  ( 1 .257  x  10"^  m^ ). 


Plunger  solenoid  actuator 


Fig.  46  Magnetic  force  versus  gap 
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Fig.  46  shows  magnetic  force  versus  gap  between  the  armature  and  the  stopper. 
The  analytical  solution  plotted  was  obtained  using  the  reluctance  method  for  the  case 
when  the  armature  is  cylindrical.  When  the  gap  is  changed  from  2mm  to  10mm,  the 
magnetic  force  decreases.  Structured  mesh  with  the  same  number  of  nodes  and 
elements  was  used  for  each  analysis  as  the  gap  was  varied.  For  the  cylindrical  plunger 
solenoid,  the  number  of  nodes  is  10984.  For  the  brick  plunger  solenoid,  the  total  number 
of  nodes  is  10876.  The  computed  values  are  higher  than  the  analytical  values  because 
the  reluctance  method  ignores  fringe  effects.  Magnetic  force  of  the  brick  plunger  is  a 
little  higher  than  the  force  of  the  cylindrical  plunger. 

2D  Clapper  solenoid  actuator  with  cantilever  beam 

A  cantilever  beam  is  attached  to  the  top  surface  of  the  armature.  The  cantilever 
beam  is  a  beam  fixed  at  one  end  and  attached  to  the  moving  armature  at  the  other  end. 
The  attachment  location  varies  as  shown  in  Fig.  47  .  The  distances  from  the  center  axis 
to  the  tip  of  the  beam  are  10mm,  20mm,  and  28mm.  The  thickness  of  the  beam  is  5mm 
and  the  length  of  the  beam  is  80mm.  The  beam  is  made  of  aluminum  whose  material 
properties  are  E  =  69  Gpa ,  v  =  0.29  and  //,  =  !. 
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Fig.  47  Planar  clapper  solenoid  actuator  with  a  cantilever  beam.  A)  d=10mm  B)  d=20mm 

and  C)  d=28mm 

For  this  analysis,  4  node  bilinear  elements  were  used.  The  total  number  of 
elements  in  the  model  was  1550.  Fig.  48  shows  the  tip  displacement  versus  ampere- 
turns.  As  N1  increases,  the  displacement  on  the  tip  increases.  When  d=20mm  or  28mm, 
the  displacements  at  the  tip  is  much  larger  than  when  d=10mm.  However,  there  is  very 
small  difference  in  tip  deflection  between  the  two  cases  where  d=20mm  and  d=28mm. 
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Fig.  48  Displacement  on  the  tip  versus.  Nl  according  to  the  location  of  the  beam 

3D  Plunger  solenoid  actuator  with  structures 

The  plunger  solenoid  actuator  with  one  of  the  three  plates  shown  in  Fig.  49 
attached  to  its  armature,  was  modeled  using  3D  elements.  The  top  views  of  each  the 
plate  structures  and  dimensions  are  shown  in  the  figure.  The  first  structure  is  a  solid 
plate,  the  second  is  a  plate  with  one  hole,  and  the  third  one  is  a  plate  with  two  holes. 
The  thicknesses  of  all  plates  are  2mm.  All  of  the  three  structures  is  attached  on  the  top 
surface  of  the  plunger  armature  as  shown  in  Fig.  50  .  The  structures  are  made  of 
aluminum. 


Round:  5mm  Round:  5mm 


Fig.  49  Top  views  of  structures  attached  on  the  piunger  armature.  A)  Soiid  piate,  B)  piate 

with  one  hole,  and  C)  plate  with  two  holes 
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Fig.  50  3D  plunger  solenoid  actuators  with  attached  structures.  A)  Solid  plate,  B)  plate 

with  one  hole,  and  C)  plate  with  two  holes 


When  NI=800  and  the  fixed  boundary  conditions  are  applied  on  the  surfaces  of  the 
structure  at  R=70mm,  the  magnetic  force  from  the  actuator  results  in  the  deflection  of 
the  plate  structures  as  shown  in  Fig.  51  . 
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Fig.  51  Deformed  structures:  A)  Solid  plate,  B)  plate  with  one  hole,  and  C)  plate  with  two 

holes 


The  maximum  displacements  for  the  structures  are  7.24x1  o  ’ mm  in  the  solid  plate, 
1.086x10“*’ mm  in  the  plate  with  one  hole,  and  8.539x10“’ mm  in  the  plate  including  two 
holes.  Fig.  52  shows  the  maximum  displacement  versus  N1  (ampere-turns)  for  each 
plate.  Among  the  three  structures,  the  plate  with  one  hole  has  the  largest  deformation 
while  the  solid  plate  has  the  highest  stiffness. 
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Fig.  52 


Maximum  displacement  versus  Nl 


7.  Discussion 

Magnetostatic  analysis  using  structured  mesh  has  been  demonstrated  using 
solution  structures  for  applying  boundary  and  interface  conditions.  The  method  allows 
geometry  imported  from  solid  modeling  systems  to  be  directly  used  for  analysis 
without  approximation  using  a  mesh.  The  approach  has  been  demonstrated  for  both  2D 
and  3D  magnetostatic  models.  Structured  mesh  is  easy  to  generate  and  the  elements  are 
regular  and  not  distorted  as  in  traditional  finite  element  mesh.  Furthermore,  the  internal 
elements  are  identical  to  each  other  and  have  the  same  stiffness  matrix  thus  reducing 
the  computation  required.  The  global  equations  were  solved  using  direct  solvers  rather 
than  iterative  solvers  whose  performance  deteriorates  when  the  global  matrix  is  ill- 
conditioned. 

Magnetic  forces  were  computed  by  integrating  the  magnetic  force  density. 
Alternative  approaches  such  as  the  virtual  work  principle  can  also  be  implemented  with 
the  implicit  boundary  approach.  If  the  force  generated  by  magnetic  forces  significantly 
deform  the  structures  on  which  the  magnets  are  mounted  then  it  is  necessary  to 
perform  a  fully  coupled  magneto-elastostatic  analysis  to  computed  the  resultant 
deformation.  In  this  case,  since  the  magnetic  field  can  be  altered  by  the  deformation  of 
the  structure,  it  leads  to  a  nonlinear  formulation  of  the  coupled  problem.  The  approach 
presented  in  this  report  needs  to  be  extended  in  the  future  to  such  nonlinear  problems 
as  well  as  dynamics  problems.  Since  3D  analysis  is  computationally  expensive,  it  is  not 
always  desirable  to  use  a  uniform  mesh  everywhere.  Mesh  refinement  techniques  are 
needed  that  locally  refines  the  mesh  while  still  using  a  structured  mesh  that  has  regular 
shaped  elements.  Development  of  such  local  refinement  techniques  is  part  of  the  future 
work  needed  to  make  3D  magnetostatic  analysis  more  feasible. 
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V.  MODELING  SHELLS  USING  IBFEM 


1.  Overview 

Shell  theories  are  extended  from  the  traditional  plate  theories.  The  two  plate 
theories  that  have  been  extensively  studied  and  implemented  in  FEM  are  the  classical  or 
Kirchoff's  plate  theory  and  the  first  order  or  Mindlin  plate  theory.  The  classical  plate 
theory  is  also  called  the  thin  plate  theory  because  it  is  based  on  Kirchhoff's  assumption 
that  for  thin  shells  shear  strain  energy  is  negligible.  Mathematically,  this  assumption 
implies  that  the  rotation  at  any  point  in  plate  is  the  slope  of  deflection  at  that  point.  The 
weak  form  of  the  plate  equations  derived  using  this  assumption  contains  volume 
integral  of  the  second  derivatives  of  deflection.  As  a  result,  the  test  and  trial  functions 
constructed  for  this  weak  form  must  be  C'  continuous  interpolation  or  approximation 
of  the  nodal  values.  This  is  difficult  to  achieve  using  traditional  interpolation  schemes 
and  therefore  the  results  obtained  using  these  elements  are  not  reliable  for  thick  plates 
and  shells.  Due  to  these  drawbacks,  the  thick  plate  theory  or  first  order  plate  theory, 
based  on  Mindlin's  assumption  [70],  are  often  preferred.  In  this  theory,  the  deflection 
and  rotation  at  any  point  in  the  plate  are  independent  of  each  other.  The  shear  strain  is 
assumed  to  be  constant  and  equal  to  the  difference  between  the  slope  and  rotation.  This 
assumption  allows  the  plate  elements  to  be  C"  continuous  because  the  weak  form 
derived  using  Mindlin's  assumption  involves  only  the  first  derivatives  of  the  deflection. 
Plate  elements  can  be  extended  to  simple  flats  shells  by  incorporating  in-plane  strains, 
allowing  it  to  handle  both  in-plane  forces  as  well  as  bending  forces.  Curved  structures 
are  approximated  by  adequate  number  of  flat  elements  [71].  The  flat  shell  elements  uses 
linear  or  bilinear  approximation  which  is  typically  C°  continuous.  More  advanced  shells 
based  on  curvilinear  coordinates  have  also  been  formulated  based  on  Mindlin  theory.. 
Straightforward  implementation  of  Mindlin  plate  theory  is  good  for  thick  plates  but 
under-predicts  deflections  for  relatively  thin  plates.  Very  thin  plates  or  shells,  exhibit 
the  so  called  shear-locking  in  plates  and  shells  [72]-[75].  Various  techniques  have  been 
proposed  to  circumvent  shear  locking.  Reduced  integration  is  the  best-known  solution 
for  shear  locking  [76].  Increasing  the  order  of  the  polynomials  used  for  the  interpolation 
also  improves  the  solution.  A  variety  of  mixed  formulations  have  also  been  successfully 
employed  to  eliminate  locking  and  to  explain  why  reduced  integration  works  [75]-[77]. 

The  analysis  of  shells  using  3D  structured  mesh  is  presented  in  this  chapter.  A 
structured  mesh  is  a  non-conforming  mesh  in  which  all  the  elements  are  regular  in 
shape  (rectangles/cuboids).  In  traditional  FEM,  shell  elements  are  two  dimensional  flat 
or  curved  elements  that  approximate  the  geometry  as  shown  in  Fig.  53  .  Fig.  53  (b) 
shows  a  shell  modeled  using  a  3D  structured  mesh  where  the  shell  geometry, 
represented  using  parametric  surfaces,  passes  through  the  elements  of  the  grid. 
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Therefore,  the  geometry  is  independent  of  the  mesh  and  could  be  represented  using  the 
exact  equations  imported  from  CAD  models.  The  use  of  structured  mesh  also  allows 
uniform  B-spline  basis  functions  to  be  used  to  represent  the  solution  so  that  the 
displacement  field  is  represented  as  C‘  or  continuous  function. 


Fig.  53  FEM  mesh  versus  structured  mesh  for  shells 

A  structured  mesh  is  much  easier  to  generate  than  a  conforming  finite  element 
mesh.  Automatic  mesh  generation  algorithms  for  FEM  can  be  unreliable  for 
complicated  geometries  often  resulting  in  poor  or  distorted  elements  and  such  distorted 
elements  is  one  of  the  main  causes  of  errors  in  FEM  solution.  Significant  amount  of  user 
intervention  is  required  rectify  the  mesh.  Moreover,  the  analysis  geometry  is  poorly 
approximated  using  elements,  especially  when  flat  shell  elements  are  used.  These 
limitations  associated  with  mesh  generation  can  be  avoided  by  using  a  structured  mesh 
for  analysis. 

The  nodes  of  a  structured  mesh  are  not  guaranteed  to  lie  on  the  analysis 
boundary,  and  therefore  the  traditional  methods  used  in  FEM  for  applying  essential 
boundary  conditions  cannot  be  used.  Implicit  boundary  method  [42]-[45]  has  been 
shown  to  be  an  efficient  and  accurate  method  for  applying  boundary  conditions  when 
the  equations  of  the  boundary  are  available.  We  extended  this  method  to  handle 
boundary  conditions  for  shells.  The  performance  and  capabilities  of  IBFEM  shell 
elements  are  studied  by  solving  several  shell  problems  in  structural  analysis. 

2.  Implicit  Boundary  Method 

The  implicit  boundary  method  uses  approximate  step  functions  of  the  boundary 
(t)(x)  =  0  to  construct  solution  structures  for  the  solution  such  that  the  displacement 
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boundary  conditions  are  guaranteed  to  be  satisfied.  Let  u-  be  displacement  component 
that  must  satisfy  the  boundary  condition  u.  =  a.  along  a  subset  of  the  boundary  of  the 
domain.  If  D  (x)  =  0  is  an  implicit  equation  of  the  boundary  r^,  then  the  solution  is 
constructed  as  [44]-[45]. 

M.(x)  =  Z),(x)w/(x)  +  i/;(x)  (5.1) 


The  solution  structure  for  the  solution  defined  in  Eq.  (5.1)  is  then  guaranteed  to 
satisfy  the  condition  m,  =  u“  along  the  boundary  defined  by  the  implicit  equation 
Z),  (x)  =  0 .  The  boundary  value  function,  u“ ,  must  be  defined  such  that  at  the  boundary 
it  has  the  prescribed  value  m"  =  a,. .  The  variable  part  of  the  solution  structure  is  the 
function  uf(x)  which  is  defined  by  piece-wise  interpolation  or  approximation  over  the 
elements  of  the  grid.  The  functions,  Z).(x),  referred  to  here  as  the  D-functions,  are 
constructed  as  approximate  step  functions,  if  essential  boundary  conditions  are 
imposed  on  the  i**^  component  of  displacement.  These  functions  are  defined  using 
implicit  equations  of  the  boundary,  ^(\)  =  0,  as 


0 


1 


^(x)  <  0 
0  <  (^(x)  <  S 
(^(x)  >  S 


(5.2) 


If  no  boundary  conditions  are  applied  on  a  boundary,  then  we  set  Z),(x)  =  1  at  that 
boundary.  The  D-functions  take  a  value  of  zero  and  has  non-zero  gradient  at 
boundaries  In  the  limit  as  this  D-function  approximates  the  Heaviside  step 

function  which  has  a  unit  value  inside  the  domain  where  ^>0  and  is  zero  at  the 
boundary  and  outside  ((z5<0).  The  use  of  an  approximate  step  functions  as  the  D- 
functions  in  implicit  boundary  method  ensures  that  only  the  boundary  elements  are 
affected  by  this  function.  The  variable  part  of  the  solution  structure  uf(x)  can  be  defined 
using  B-spline  or  other  approximations  that  provide  a  high  order  continuity.  Uniform 
B-spline  elements  applicable  to  structured  mesh  have  been  studied  in  the  past  [44],  [46] 
to  demonstrate  that  the  implicit  boundary  method  can  be  used  to  apply  boundary 
conditions  even  when  the  approximation  used  does  not  have  Kronecker's  delta 
properties.  In  other  words,  even  though  the  nodal  values  are  not  equal  to  the  B-spline 
approximation  values  at  the  nodes,  the  solution  structure  in  Eq.  (5.1)  guarantees  that  the 
essential  boundary  conditions  are  imposed. 

Eor  shells,  the  boundary  is  a  curve  on  a  surface  and  essential  boundary  conditions 
specify  values  of  displacements  along  the  curved  boundary  or  the  local  slopes  on  these 
curves.  Both  these  boundary  conditions  can  be  imposed  by  appropriately  constructing 
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the  characteristic  functions  4>{\)  for  the  boundary  curves.  To  apply  displacement 
boundary  conditions,  where  the  slope  is  not  fixed,  the  characteristic  function  is  defined 
as  the  radial  distance  from  the  curve.  In  order  to  model  clamped  boundary  condition, 
where  both  slope  and  displacement  are  fixed,  the  characteristic  function  is  defined  as 
the  distance  from  a  imaginary  surface  that  passes  through  the  boundary  curve  and  is 
normal  to  the  shell  surface. 

The  weak  form  of  the  elastostatic  boundary  value  problem  is  the  principle  of 
virtual  work  which  can  be  written  in  the  following  form: 

=  |8u^tfi?r  + (5.3) 

Q  Q 

Using  the  solution  structure  in  Eq.  (5.1)  the  stresses  and  strains  can  also  be 
decomposed  into  a  homogeneous  part  and  a  boundary  value  part. 

£  =  +  £'' 

(5.4) 

0  =  C£  =  C(£* +£'')  =  «*+ o’* 


Where, 
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The  Dirichlet  function,  D. ,  is  defined  as  per  Eq.  (5.2)  for  each  displacement 
component  for  which  essential  boundary  condition  is  specified.  Substituting  the 
preceding  equations  into  the  weak  form,  the  following  relation  is  obtained,  where  the 
known  quantities  have  been  moved  to  the  right  hand  side  of  the  equation. 

I  (5.5) 

Q  Q  Tf  Q 

The  boundary  value  function  u-  must  be  constructed  such  that  at  the  boundary 
where  Dirichlet  boundary  conditions  are  applied  it  must  have  value  equal  to  the 
specified  boundary  condition  at  that  boundary.  When  Dirichlet  boundary  conditions 
are  specified  at  multiple  boundaries  and  the  assigned  values  are  not  the  same  for  all 
boundaries,  the  boundary  value  function  must  take  on  appropriate  values  at  respective 
boundaries  and  transition  smoothly  in  between.  To  accomplish  this,  the  boundary  value 
function  is  constructed  by  piecewise  interpolation  or  approximation  using  the  shape 
functions  of  the  elements  through  which  the  boundary  is  passing.  A  displacement 
component  u“  is  defined  within  an  element  as  where  are  the  nodal 

J 

values  of  the  displacement  component. 
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3.  B-Spline  approximations 

B-spline  shape  (basis)  functions  are  traditionally  constructed  using  a  recursive 
definition  (Farin,  2002).  The  parameter  space  is  partitioned  into  elements  using  a  knot 
vector  (equivalent  to  a  collection  of  nodes).  General  methods  are  available  to  insert 
knots  and  elevate  the  order  of  the  polynomial.  However,  for  implementation  in 
structured  mesh  finite  element  analysis,  it  is  convenient  to  define  shape  functions 
within  each  element  of  the  mesh  such  that  the  parametric  domain  of  the  element  is  [-1,1] 
in  the  three  parametric  directions.  The  use  of  such  shape  functions  for  such  uniform  B- 
splines  shape  functions  have  been  demonstrated  in  past  work  for  3D  elastostatic 
elements  [44].  Here  we  briefly  summarize  the  main  ideas  and  provide  the  shape 
functions  that  were  used  for  shell  analysis. 

The  shape  functions  for  one  dimensional  elements  are  derived  using  the 
continuity  requirements  between  neighboring  elements.  Any  k‘^  order  B-spline  has 
A: +  1  support  nodes.  A  quadratic  B-spline  element  in  one-dimension  has  three  nodes  and 
therefore  it  is  represented  by  three  shape  functions.  The  coefficients  of  the  shape 
function  are  determined  by  requiring  the  approximations  and  its  slope  to  be  continuous 
between  elements.  Additionally,  the  shape  functions  are  required  to  be  partition  of 
unity.  The  expression  of  the  shape  functions  are  given  below, 

N,=\{\-2r-r^) 

O 

7V,=i(6-2r^)  (5.6) 

O 

Afj  =-(l  +  2r  +  r") 

8 

The  plots  of  the  shape  functions  in  Eq.  (5.6)  are  shown  in  Fig.  54 


Fig.  54  Shape  function  of  one  dimensionai 
quadratic  B-spiine  eiement 

A  cubic  B-spline  element  in  one-dimension  has  four  nodes  and  therefore  it  is 
represented  by  tour  shape  functions.  These  shape  functions  are  derived  by  requiring  the 
approximation,  its  slope  and  curvature  to  be  continuous  between  elements.  These  shape 
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functions  are  also  required  to  be  a  partition  ot  unity.  The  expression  ot  the  shape 
functions  are  given  below, 

N,  =— (l-3r  +  3r"-rh 
‘48 

N,  =— (23-15r-3r"+3rh 

48 

(5.7) 

N,  =— (23  +  15r-3r"-3rh 
^  48 

Af,  =— (l  +  3r+3r"+rh 

^  /I  O  ^  ^ 


The  plots  of  the  shape  functions  in  Eq.  (5.7)  are  shown  in  Fig.  55  This  figure  shows 
that  the  B-splines  are  not  unity  at  respective  nodes  and  do  not  vanish  at  other  nodes. 
This  shows  that  they  do  not  interpolate  nodal  values.  Instead,  smooth  approximations 
of  the  nodal  values  are  created. 


Fig.  55  Shape  function  of  one  dimensional  cubic  B-spline  element 

The  shape  functions  for  the  higher  dimensional  B-spline  elements  are  constructed 
as  products  ot  the  shape  functions  tor  one-dimensional  B-splines.  The  structured  grid 
consists  of  regular  hexahedra  (cube/cuboid),  so  the  mapping  from  parametric  space  to 
the  physical  space  is  linear.  The  mapping  between  parametric  and  physical  space  can  be 
defined  as 

(5.8) 

'  2  '  2  ' 

where,  x'  and  x“  are  the  lower  and  upper  bounds  respectively  for  nodal 
coordinates  of  any  given  element. 

The  shape  functions  for  two  and  three-dimensional  cubic  B-splines  are  constructed 
as  a  product  ot  one-dimensional  corresponding  one-dimensional  shape  functions.  The 
shape  functions  can  be  expressed  as, 

=  ^,.(r)iV/5)7V,(0 
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4.  Numerical  Implementation 

The  B-spline  shell  elements  are  implemented  as  three-dimensional  regular 
elements  (cube/cuboid)  that  are  part  of  a  structured  mesh  and  each  element  contains  a 
portion  of  the  surface  representing  the  shell  as  shown  in  Fig.  53  The  stiffness  matrix  for 
linear  elastic  elements  is  computed  as 

K‘  =  \[BY[C][BW  (5.10) 


Only  a  portion  of  the  shell  surface  that  intersects  with  the  element  contributes  to 
its  stiffness.  The  volume  integration  in  Eq.  (5.10)  is  evaluated  as  a  combination  of  area 
and  thickness  integration.  In  order  to  integrate  accurately  over  the  portion  of  the  surface 
passing  through  the  element,  intersection  between  the  surface  and  the  element  is 
computed  and  triangulated.  The  stiffness  is  then  computed  by  integrating  over  these 
triangles  and  adding  the  results  as  follows 


nt  2 


[K]  =  Zj  j[sf[C][SM4 


m 


(5.11) 


In  Eq.  (5.11)  4  is  the  area  of  the  i*  triangle  and  h  is  the  thickness  of  the  shell. 

For  boundary  elements  that  contain  boundaries  with  essential  boundary 
conditions,  the  D-functions  sharply  dip  from  unit  value  to  zero  near  the  boundary  over 
a  narrow  band  of  width  The  gradients  of  D-function  over  these  narrow  bands 

can  have  large  values.  In  order  to  accurately  compute  the  contribution  of  the  D-function 
and  its  gradient  to  the  stiffness  matrix,  strain-displacement  matrix,  [B]  is  split  into  two 
parts  as  [B'\  =  {B^\+{B^'\,  such  that  [5,]  contains  only  derivatives  of  the  shape  functions 
while  contains  derivatives  of  the  D-function  [42]-[44].  Away  from  the  boundaries 
with  specified  boundary  conditions,  Z),  is  unity  and  therefore  [^4  =  0  because  the 
derivative  of  Z),  are  zero.  Now  the  stiffness  matrix  can  be  expressed  as 

A 

[ZZ],  =  ]  [C]{{B,-\  +  [B,-\)dA,dh  (5.12) 

_*  4- 


Expanding  this  expression,  it  can  be  seen  that  the  stiffness  matrix  has  the 
following  three  components. 


=[z:,]+([z:4+[z:,f)+[z:3] 


(5.13) 
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[^,]=f  \[B,f [C][B,YIA 


V 


Vih 


[K,]  =  f  j [B,f [C][B,^flA,  +  j [S,f [C][5,M4 


V 


to 


[/:3]=J  \[B,f [C][b,yia 


to 


(5.14) 


(5.15) 


(5.16) 


[.s:,]  must  be  computed  by  sub-dividing  the  surface  into  triangles.  and  [.STj] 
contain  terms  involving  gradients  of  D-functions  and  zero  everywhere  except  in  the 
vicinity  of  boundaries  that  have  essential  boundary  conditions  specified.  Using  this  fact 
the  computation  of  [K^]  and  [a:,]  can  be  converted  into  line  integrals  along  the  boundary 
curves  as  shown  below. 


[^2]  =  ZJ  JJtoftoll^J^  dcjxiljdh 

j=^  _hlj  <!> 


(5.17) 


In  Eq.  (5.17),  n^  is  the  total  number  of  line  segments  /.  used  to  represent  the 
boundary  curve  within  the  element  e .  Similarly,  the  last  term  can  be  written  as 


\.K,-\  =  Y^]\\{BJ[C][B,-\  ■^-^^d(/)dljdh 


7=1  _h  I  ^ 


(5.18) 


The  boundary  load  terms  are  computed  again  by  approximating  the  shell 
boundaries  within  each  element  as  a  set  of  line  segments  and  evaluating  the  load 
integral  along  each  line  using  Gauss  quadrature  as 

h 

{/}  =  ill  {NY  {t}dTdh  (5.19) 

i=\  ^  1. 

2 

In  Eq.  (5.19),  {t}  is  the  traction  acting  on  the  boundary,  /,  are  the  line  segments  and 
the  shell  thickness  is  h .  Pressure  loads  and  body  forces  require  integration  of  the  shell 
surface  which  must  be  approximated  as  triangles  for  integration  purpose.  For  example, 
body  force  contribution  is  evaluated  as 
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{F‘‘}=^^{NY{b}dAdh 


(5.20) 


-h  A 


2 


In  Eq.(5.20),  the  body  forces  are  evaluated  as  a  combination  of  through  thickness 
and  area  integrals.  The  integration  over  area  is  computed  as  a  summation  of  the 
integration  over  triangles  generated  on  the  surface. 


5.  Results  and  Discussion 


In  this  section,  the  3D  structures  mesh  approach  to  shell  analysis  is  demonstrated 
using  several  examples  and  the  results  obtained  by  this  approach  has  been  referred  to 
as  IBFEM  results.  Quadratic  B-splines  have  27  nodes  are  referred  to  as  IBFEM27N  and 
cubic  B-spline  elements  are  IBFEM64N. 

Centrally  Loaded  Square  Clamped  Plate  in  Bending 

A  square  plate  as  shown  in  Fig.  56  is  loaded  by  a  uniform  pressure  of  100  psi.  The 
material  used  is  structural  steel.  The  problem  description  including  dimensions  and 
loads  are  shown  in  the  figure.  The  plate  is  clamped  along  all  its  edges.  Both  the  64-node 
cubic  B-spline  elements  and  27-node  quadratic  B-spline  elements  were  used  for  the 
analysis  and  compare  with  results  obtained  using  traditional  shell  elements 
implemented  in  commercial  software  (ABAQUS). 


L  b  =  10  in 
r  =  0.1  in 

F  =  100  psi  (downward  pressure) 

Young's  Modulus  =  1.1 1966  X  10‘®  psi 
Poison's  ratio  =  0.27 


Fig.  56  Centrally  loaded  square  plate 
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Fig.  57  IBFEM  results  (using  cubic  B-spline  64  node  elements  (5x5  Mesh) 


77 


Fig.  58  Abaqus  result  using  S4R  (10x10  Mesh) 

Table  IV.  Square  plate:  results  for  vertical  displacement  at  the  middle  of  the 
square  plate,  based  on  various  meshes  and  element  types 


Element  type 

Mesh 

Max.Vertical 

displacement 

X  10'^  (in) 

IBFEM  64N 

2x2 

3.562 

5x5 

4.051 

10  X  10 

4.055 

IBFEM  27N 

lOx  10 

3.966 

20  X  20 

4.043 

40  x40 

4.059 

ABAQUS  S4R 

5x5 

3.703 

10  X  10 

4.054 

20x20 

4.055 

The  results  obtained  using  various  elements  are  listed  in  table  IV.  Using 
traditional  FEM  shell  elements  (ABAQUS  S4R  elements),  the  maximum  deflection  at  the 
center  numerically  converges  to  a  value  of  4.055  x  10"^  inch  of  vertical  displacement  at 
the  middle  of  the  square  plate.  The  table  shows  that  the  results  of  IBFEM  64  node  shell 
elements  converge  to  the  same  answer  with  fewer  number  of  elements.  These  elements 
use  cubic  polynomial  representation  of  the  solution  and  hence  represent  shear  strains 
through  thickness  as  parabolic  and  hence  are  able  to  represent  the  exact  solution  more 
closely.  The  quadratic  27  node  B-spline  shell  elements  require  a  much  larger  number  of 
elements  to  reach  the  same  solution. 

Micro  Air  Vehicle  Wing  design 

The  dimensions  of  the  structure  for  this  problem  are  shown  in  Fig.  59  .  The  wing 
is  subjected  to  a  normal  pressure.  One  end  is  fixed  to  the  fuselage  while  the  side  edges 
are  free.  The  material  properties,  loads  and  constraints  are  shown  in  the  figure.  The 
thickness  of  the  shell  is  0.1  inch.  Results  computed  using  64  node  B-spline  elements  and 
traditional  FEM  (ABAQUS  S8R5)  are  compared  here. 
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Yoiin8'&  Modulus  =  3.0  X 
Poison's  ratio  =  0.3 


Fig.  59  Micro  air  vehicle  wing 
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Fig.  60  Displacement  plot  (IBFEM  64N)  for  micro  air  vehicle  wing  (10x5  Mesh) 


Fig.  61  Displacement  plot  for  micro  air  vehicle  wing  (ABAQUS  S8R5,  542  elements) 

The  figures  show  that  the  results  of  IBFEM  shell  elements  the  ABAQUS  8  node 
iso-parametric  reduced  integration  elements  are  very  close. 

Design  of  vibrating  /  osciilating  wings  for  MAV 

To  design  oscillating  structures  such  as  flapping  wings  it  necessary  to  analyze  the 
dynamics  of  structures.  Modal  analysis  computes  the  natural  frequencies  and  the  mode 
shapes  of  vibration  of  a  structure.  New  concepts  for  oscillating  wing  design  were 
studied  using  modal  analysis  capability  that  was  recently  added  to  IBFEM  software. 

Eig.  62  shows  a  structural  mechanism  that  can  oscillate  using  internally  applied 
forces  generated  using  embedded  electromagnetic  actuators.  The  analysis  was 
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performed  using  novel  shell  elements  implemented  in  IBFEM  software  that  can  perform 
the  analysis  for  relatively  complex  geometry  without  needing  a  conforming  surface 
mesh.  Fig.  62  (a)  and  (c)  show  the  3D  grid  and  the  applied  forces  that  were  used  for  the 
analysis.  The  deformed  shape  of  the  wing  like  geometry  is  shown  in  Fig.  62. 


Fig.  62  Oscillating  structural  mechanism 

The  modes  of  vibration  of  this  structure  were  computed  to  understand  the  shapes 
in  which  this  structure  can  be  made  to  vibrate  or  oscillate  by  inducing  resonance.  Fig.  63 
shows  the  four  modes  of  vibration  of  this  structure. 


80 


(c)  Mode-Ill  (d)Mode-IV 


Fig.  63  Mode  shapes  of  vibration 

Modes  1-,  II  and  IV  are  similar  to  cantilever  beam  oscillations  but  Mode-Ill 
involves  twisting.  Clearly  none  of  these  modes  are  ideally  suited  for  producing  lift  and 
thrust  for  flapping  flight  for  MAVs.  Therefore,  the  design  challenge  is  to  compute  the 
right  shape  and  reinforcement  for  the  wing  such  that  the  first  mode  of  vibration  will 
produce  a  fanning  action  that  produces  thrust.  This  will  be  attempted  in  future  designs 
using  intuitive  changes  to  the  design.  In  future  work,  it  is  hoped  that  we  will  develop 
topology  optimization  techniques  to  compute  appropriate  reinforcements  for  a  desired 
mode  of  vibration. 


6.  Discussion 

Shell  elements  based  on  uniform  B-spline  shape  function  were  studied  and  shown 
to  be  accurate  when  compared  with  traditional  FEM  shell  elements  implemented  in 
commercial  software.  One  of  the  key  advantages  of  using  these  elements  is  that  a 
structured  mesh,  which  is  easy  generate  automatically,  can  be  used  for  the  analysis. 
Moreover,  the  geometry  is  represented  accurately  using  equations  rather  than 
approximated  by  the  mesh.  Implicit  boundary  method  for  applying  boundary 
conditions  is  extended  to  shells.  Both  quadratic  and  cubic  B-spline  shape  functions  were 
tested  and  it  was  found  that  cubic  elements  provided  very  good  results  with  fewer 
elements  than  quadratic  elements.  Computational  cost  is  higher  for  these  elements 
compared  to  traditional  shell  elements  but  often  fewer  elements  are  needed  to  get 
accurate  results  with  cubic  elements.  The  time  taken  to  create  the  model  is  significantly 
lower  because  structured  mesh  generation  is  easily  automated.  Future  work  in  this  area 
includes  extending  these  elements  to  non-linear  applications  involving  large  elastic 
deformation  and  elasto-plastic  deformations.  The  performance  of  the  IBFEM  shell 
elements  can  be  further  improved  by  incorporating  adaptive  mesh  refinement.  The 
stress  concentrations  problems  can  be  analyzed  faster  with  adaptive  mesh  refinement. 
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The  shell  elements  using  IBFEM  can  also  be  extended  to  mixed  solid-shell  type  analysis 
where  thin  shells  are  attached  to  solid  structures. 
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VI.  DESIGN  OF  ACTUATORS  FOR  FLAPPING  WING 


1.  Design  criteria 

In  this  chapter,  IBFEM  is  used  as  an  actuator  design  tool.  Solenoid  actuators  with 
clapper  armature  or  plunger  armature  and  coil  actuators  are  designed  under  given 
design  criteria.  The  main  application  of  interest  is  actuators  for  flapping  wing  micro-  air 
vehicles.  The  micro  air  vehicles  are  small  size  aerial  vehicles  with  a  wing  span  of  about 
15  cm.  For  a  micro  air  vehicle  actuator,  three  main  design  criteria  were  used:  coil 
winding  area,  iron  weight  and  size. 

To  illustrate  the  analysis  and  design  methodology  using  IBFEM,  the  process  is 
demonstrated  here  using  a  set  of  assumed  specification.  An  acceptable  value  for  coil 
winding  area  was  assumed  to  be  I2mm^ .  The  entire  coil  winding  area  cannot  be 
assumed  to  be  only  copper  because  the  coil  includes  copper  and  other  materials 
including  insulator  and  air  gaps.  The  relation  between  the  coil  winding  area  and  the 
copper  conductor  is  defined  as  the  packing  factor  [68].  When  the  coil  is  wound 

tightly,  the  packing  factor  can  be  up  to  75% .  Thus,  the  packing  factor  is  assumed  to  be 
70% .  The  area  of  the  coil  wire  is  =  5.01xl0“^ww^  when  the  number  of  the  American 
Wire  Gauge  (AWG)  is  40.  Therefore,  the  number  of  coil  turns  N  is  given  as  following 
equation 

F  S 

A  =  ^^  =  1680  (6.1) 

When  the  coil  is  wound  to  a  cylindrical  bobbin  with  the  diameter  of  4  mm,  the  total 
length  of  the  coil  becomes  about  21.1  m.  As  the  resistance  per  meter  for  the  given  AWG 
is  3.44  Q/w,  the  total  resistance  of  the  coil  is  approximately  72  Q.  In  this  research,  it  is 
assumed  that  the  coil  wire  can  allow  current  of  100mA  to  flow.  The  allowable  current  of 
the  wire  is  based  on  plastic  insulation.  In  case  of  the  current  source,  LT3092, 
manufactured  by  Linear  technology,  the  maximum  output  current  is  200mA.  The  input 
voltage  range  is  from  1.2  V  to  40  V  so  that  LT3092  can  be  operated  by  using  a  small  size 
battery.  Thus,  the  current  of  100mA  can  be  obtained  using  LT3092.  If  N=1680  and  the 
amount  of  current  1  can  be  controlled  by  a  digital  processor,  N1  (ampere-turns)  can  vary 
from  0  to  168 .  N1  is  also  called  magnetomotive  force.  In  order  to  compare  several 
designed  actuators  in  the  later  section,  the  magnetomotive  force,  NI=30,  is  used.  When 
the  maximum  current  (100mA)  flows  in  the  coils,  the  coils  can  have  maximum  energy 
dissipation  as  W  =  I^R  =  0.72  [W] . 
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Another  design  criterion  is  the  iron  weight  of  the  actuator.  The  iron  core 
contributes  most  to  the  weight  of  the  actuator.  Iron  has  high  relative  permeability  so 
that  using  iron  for  the  core  of  the  actuator  can  intensify  the  magnetic  force  of  the 
actuator.  Therefore,  there  is  a  trade-off  between  the  iron  weight  of  the  actuator  and  the 
usage  of  the  iron.  As  the  second  design  criterion,  we  assumes  that  the  total  iron  weight 
can  be  up  to  lOg.  In  order  to  compute  the  iron  weight  of  a  actuator  design,  the  mass 
density  of  iron  is  assumed  to  be  7874A:g/w^  =  7.874 xlO”^g/mw^  .  The  third  design 
criterion  is  the  size  of  the  actuator.  In  order  to  install  an  actuator  inside  the  micro  air 
vehicle,  we  have  assumed  that  the  actuator  should  fit  inside  of  a  cube  with  edge  length 
of  2cm. 

Here  we  first  compare  several  actuator  designs  by  computing  the  force  generation 
capacity  using  the  magneto-elastostatic  analysis  capability  implemented  in  IBFEM.  Both 
solenoid  type  actuators  and  coil  actuators  were  studied. 

2.  Solenoid  actuators  with  ciapper  armature 

A  solenoid  actuator  with  a  clapper  armature  is  shown  in  Fig.  64  .  The  clapper 
armature  can  move  in  the  y-direction  when  the  coils  carry  current.  The  armature  and 
the  stator  are  made  of  thin  steel  laminations  with  relative  permeability  of  2000,  a  thin 
steel  lamination  that  can  reduce  heating  caused  by  eddy  current.  When  the  current 
flows  through  a  coil,  magnetic  flux  is  created  as  shown  in  Fig.  64  .  The  magnetic  flux 
follows  a  closed  path.  If  magnitude  of  the  magnetic  flux  is  known,  a  magnetic  force  can 
be  estimated  using  the  Maxwell  stress  tensor  method. 


Fig.  64  Solenoid  actuator  with  clapper  armature 

We  assume  the  dimensions  are  /I  =  /3  =  0.5mm  ,  12  =  15  =  3mm  ,  14  =  16  =  6.5mm  , 
g  =  Imm ,  and  wl  =  w2  =  w3  =  1mm  .Using  the  reluctance  method,  the  magnetic  in  the  air 
gap  is  estimated  to  be  approximately  (f)  =  6.285  x  10“’  x  M[Wb] .  Using  this  estimate,  the 
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magnetic  flux  density  and  the  magnetic  field  density  is  =6.285x10  ''xM[T]  and 
^,a,=500xM[A/m] 

Using  Maxwell's  stress  tensor  method,  the  normal  magnetic  pressure  can  be 
obtained  as 


F  =  f  (4 X  wl)  =  6.283 x  10'"  x (M)"  [N] 


(6.2) 


When  M  =  30,  the  magnetic  force  per  unit  length  becomes  0.565  N/m.  The  magnetic 
flux  and  the  magnetic  field  density  are  1.886x10“^  [T]  and  1.5x10"*  [A/m] .  Using  the 

same  geometry  of  the  actuator.  Fig.  65  (A)  shows  contour  plot  of  magnetic  vector 
potential  computed  using  IBFEM.  Fig.  65  (B)  shows  magnetic  field  density  in  the  y- 
direction.  The  computed  magnetic  field  density  in  the  gap  is  approximately  equal  to  the 
analytical  value  of  1.5x10"*  [A/m] .  For  this  analysis,  the  asymptotic  boundary  conditions 


are  applied  on  the  outer  boundaries.  The  computed  force  is  0.633  N/m.  The  computed 
force  and  the  analytical  solution  are  quite  close.  The  difference  between  them  results 
from  a  fringing  flux,  which  the  reluctance  method  ignores.  As  the  iron  domain  is  40mm^ , 
the  iron  weight  per  unit  depth  is  0.3 1 5g  /  mm  . 
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Fig.  65  Computed  results  using  IBFEM.  A)  Contour  plot  of  magnetic  vector  potential,  and 

B)  Magnetic  field  density 

Based  on  the  initial  design  shown  in  Fig.  64  ,  three  different  designs  are  created  by 
changing  the  geometry  of  the  armature  and  the  stator.  Fig.  66  shows  three  designs  with 
dimensions.  The  dimensions  shown  are  all  in  mm.  Even  though  the  geometry  is 
different  for  each  actuator,  same  mesh  densities  was  used  in  the  structured  mesh 
created  for  all  the  three  designs.  Nine  node  B-spline  elements  were  used  for  the  analysis. 
The  total  number  of  nodes  was  2586.  The  asymptotic  boundary  conditions  were  applied 
on  the  outer  boundaries.  The  magnetomotive  force,  NI  =  30,  was  applied. 
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Fig.  66  Clapper  solenoid  actuators.  A)  Design  1,  B)  Design  2,  and  C)  Design  3 

Fig.  67  shows  contour  plots  of  magnetic  vector  potentials  for  all  designs.  Based  on 
those  computations.  Table  V  shows  comparison  between  the  three  designs  in  terms  of 
the  magnetic  force  and  the  iron  weight. 
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Fig.  67  Contour  plots  of  magnetic  vector  potential  for  clapper  solenoid  actuators.  A) 

Design  1,  B)  Design  2,  and  C)  Design  3 


Table  V.  Comparison  for  three  clapper  solenoid  actuators  (Nl=30) 


Design  1 

Design  2 

Design  3 

Force  per  m 

0.691  N/m 

0.509  N/m 

0.681  N/mm 

Iron  weight  per  mm 

0.3l5g  /  mm 

0.2S3g  /  mm 

0.320g/  mm 
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According  to  Table  V,  the  first  actuator  design  produces  the  largest  force  among 
them.  The  second  actuator  is  the  lightest  one. 

3.  Solenoid  actuators  with  plunger  armature 

A  solenoid  actuator  with  a  plunger  armature  is  shown  in  Fig.  68  .  The  shape  of  the 
plunger  armature  can  be  a  brick,  cylinder,  or  conics.  There  are  two  gaps  of  g  and  gs  that 
the  flux  can  enter  and  leave  the  plunger  armature.  Useful  magnetic  force  is  produced 
only  at  g. 

Plunger  Armalute 
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Fig.  68  Solenoid  actuator  with  plunger  armature 

The  armature  and  the  stator  are  made  of  steel  with  the  relative  permeability  of  2000. 
The  dimensions  are  11  =  7 mm  ,  12  =  2.5mm  ,  13  =  0.9mm  ,  14  =  4.5mm  ,  15  =  1. 5mm  , 
16  =  3.5mm  ,  gs  =  0.lmm,  g  =  \mm  and  w\  =  w2  =  w3  =  w4  =  \mm  .  Using  the  reluctance 
method,  the  magnetic  flux  is  =  1.132x10"®  x  AT  [Wb].  Based  on  the  magnetic  flux,  the 
magnetic  flux  density  and  the  magnetic  field  density  in  each  air  gap  can  be  computed  as 

follow  =1.132x10“^  xAT[T]  and  =  900.8  x  AT  [A/m].  Using  the  magnetic 

Ao 

field  in  the  gap,  the  useful  magnetic  force  can  be  calculated  as  follow.  When  AT  =  30, 
the  force  becomes  0.918  N.  The  magnetic  flux  and  the  magnetic  field  density  are 
3.4x10"^  [T]  and  2.702 x  10"^  [A/m] .  Fig.  69  shows  the  computed  results:  contour  plot  of 

magnetic  vector  potential  and  magnetic  field  density  in  the  y-direction.  For  the  analysis, 
the  same  mesh  density  was  used  in  the  clapper  solenoid  actuator  model.  Nine  node  B- 
spline  elements  were  used.  Asymptotic  boundary  conditions  were  applied  on  the  outer 
boundaries.  The  computed  total  force  on  the  plunger  armature  is  0.9130  N/m.  As  the 
total  area  of  the  iron  is  3S.8mm^ ,  the  iron  weight  per  unit  depth  is  0.306g  /  mm  . 
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Fig.  69  Computed  results  using  IBFEM.  A)  Contour  plot  of  magnetic  vector  potential,  and 

B)  Magnetic  field  in  the  y-direction 


Modifying  the  initial  design  in  shown  in  ,  three  different  designs  are  created  as 
shown  in  Fig.  70  .  The  dimensions  shown  are  all  in  mm.  For  all  designs,  the  area  of  the 
2D  model  is  as  same  as  one  of  the  clapper  solenoid  actuator.  Nine  node  B-spline 
elements  were  used  tor  the  analysis.  The  total  number  of  nodes  was  2548.  The 
asymptotic  boundary  conditions  were  applied  on  the  outer  boundaries.  The 
magnetomotive  force,  AT  =  30,  was  used. 


Fig.  70  Plunger  solenoid  actuators.  A)  Design  1,  B)  Design  2,  and  C)  Design  3 

Fig.  71  shows  contour  plots  of  magnetic  vector  potentials  for  all  designs.  Based  on 
those  results.  Table  VI  shows  comparison  for  three  designs  in  terms  of  the  magnetic 
force  and  the  iron  weight. 
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Fig.  71  Contour  plots  of  magnetic  vector  potential  for  plunger  solenoid  actuators.  A) 

Design  1,  B)  Design  2,  and  C)  Design  3. 


Table  VI.  Comparison  for  three  plunger  solenoid  actuators  (Nl=30) 


Design  1 

Design  2 

Design  3 

Force  per  m 

1.154  N/m 

0.262  N/m 

2.801  N/mm 

Iron  weight  per  mm 

0.306g/  mm 

0.294g  /  mm 

0.3  lOg/  mm 
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According  to  Table  VI,  the  third  plunger  solenoid  actuator  model  can  produce  the 
largest  force  among  them.  The  second  actuator  model  is  the  lightest  one. 

4.  Solenoid  actuators  with  a  combined  plunger  &  clapper  armature 

A  combined  plunger  &  clapper  armature  is  an  armature  that  includes  a  clapper  & 
plunger  to  increase  the  force.  Three  designs  were  studied.  The  sizes  of  those  actuators 
are  similar  to  one  of  the  clapper  solenoid  actuator  in  the  previous  section.  Magnetic 
force  and  iron  weight  for  each  actuator  are  characterized. 

Fig.  72  shows  three  solenoid  actuators  with  a  mixture  armature.  The  dimensions 
shown  are  all  in  mm.  For  all  designs,  the  area  of  the  2D  model  is  as  same  as  one  of  the 
clapper  solenoid  actuator.  For  the  analysis,  the  same  mesh  density  is  applied  for  all  the 
designs.  Nine  node  B-spline  elements  were  used  for  the  analysis.  The  total  number  of 
nodes  was  2674.  The  asymptotic  boundary  conditions  were  applied  on  the  outer 
boundaries.  The  magnetomotive  force  NI  is  equal  to  30. 


Fig.  72  Three  combined  piunger  &  ciapper  soienoid  actuators.  A)  Design  1,  B)  Design  2, 

and  C)  Design  3 

Contour  plots  of  magnetic  vector  potentials  are  show  in  Fig.  73  .  Based  on  those 
results.  Table  VII  shows  comparison  for  three  designs  in  terms  of  the  magnetic  force 
and  the  iron  weight. 
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Fig.  73  Contour  plots  of  magnetic  vector  potential  for  combined  plunger  &  clapper 
solenoid  actuators.  A)  Design  1,  B)  Design  2,  and  C)  Design  3. 

Table  VII.  Comparison  for  three  combined  plunger  &  clapper  solenoid  actuators 
(A//=30) 

_ Design  1 _ Design  2 _ Design  3 _ 

Magnetic  force  per  0.961  AZ/at?  0A^3  N/m  1.30  AZ/m 
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0.339g/  mm 


0.340g/  mm 


m 

Iron  weight  per  mm  0.342g/mm 


According  to  Table  VII,  the  third  one  can  produce  the  largest  force  among  them.  The 
second  actuator  is  the  lightest  one. 


5.  Coil  actuators 

In  this  section,  three  coil  actuators  are  studied  under  the  given  design  criteria.  The 
coil  actuators  include  permanent  magnets  so  that  higher  magnetic  flux  can  be  produced 
without  current  flowing  though  the  actuator  resulting  in  heating  losses.  The  coil 
actuator  has  magnetic  force  in  the  conductive  coil  instead  of  one  on  the  armature  of  the 
solenoid  actuator,  a  magnetic  force  that  is  computed  using  the  Lorentz  force  equation. 

Fig.  74  shows  three  coil  actuators.  The  dimensions  shown  are  all  in  mm.  For  all  the 
designs,  the  rectangular  area  of  the  2D  model  including  the  air  domain  is  18xl8mm^. 
The  coil  winding  area  is  the  same  as  in  previous  solenoid  actuators.  Two  permanent 
magnets  have  different  direction  with  the  same  remanant  flux  of  0.8  T.  The  first  design 
is  based  on  a  typical  voice  coil  actuator  in  a  loudspeaker.  The  second  is  a  design  to 
remove  the  iron  laminates  from  the  first  design.  The  third  design  has  a  gap  at  the 
bottom  portion  of  the  iron  laminate,  which  allows  the  coil  to  move  freely  in  the 
downward  direction.  For  the  analysis,  the  same  mesh  density  is  applied  for  all  designs. 
Four  node  bilinear  elements  were  used.  The  total  number  of  nodes  was  12036.  The 
asymptotic  boundary  conditions  were  applied  on  the  outer  boundaries.  NI  is  equal  to  30. 


Fig.  74  Three  coil  actuators.  A)  Design  1,  B)  Design  2,  and  C)  Design  3 

The  magnitudes  of  the  magnetic  flux  densities  are  show  in  Fig.  75  .  As  the  Lorentz 
force  is  proportional  to  the  magnetic  flux  density,  a  stronger  magnetic  flux  density  near 
the  moving  coils  can  create  a  larger  force.  The  first  design  has  the  strongest  magnetic 
flux  near  the  moving  coils.  Based  on  those  results.  Table  VIII  shows  comparison  for 
three  designs  in  terms  of  the  Lorentz  force  on  the  moving  coil  and  the  iron  weight  of  the 
actuator. 
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Fig.  75  Magnitude  of  B  fields  for  three  coil  actuators.  A)  Design  1,  B)  Design  2,  and  C) 

Design  3. 


Table  VIII.  Comparison  for  coil  actuators  {NI=Z0) 


Design  1 

Design  2 

Design  3 

Force  per  m 

21.0  A//m 

2.64  N/m 

9.12  A//m 

Iron  weight  per  mm 

0.428g/  mm 

0  g!  mm 

0.394g/  mm 

According  to  Table  VIII,  the  first  coil  actuator  can  produce  the  largest  force  among 
them.  The  second  design  is  the  best  design  if  the  actuator  weight  is  the  most  critical 
criterion. 


6.  The  best  actuator  among  the  designed  actuators 

Using  IBFEM,  four  types  of  actuators  were  examined  such  as  the  clapper  solenoid 
actuator,  the  plunger  solenoid  actuator,  the  combined  plunger  &  clapper  solenoid 
actuator,  and  the  coil  actuator.  Among  them,  the  coil  actuator  is  the  best  actuator 
considering  the  given  design  criteria.  Among  several  designs  for  the  coil  actuators,  the 
first  design  is  used  for  flapping  wings  of  the  micro  air  vehicle  as  shown  in  Fig.  76  Coil 
actuators  have  several  advantages  over  the  solenoid  actuators  that  were  studied.  The 
most  important  being  the  force  direction  can  be  reversed  by  changing  the  direction  of 
current.  The  magnitude  does  not  vary  due  to  the  deformation  or  movement  of  the  coil 
as  long  as  the  coil  is  within  the  uniform  field  created  by  the  magnets. 
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Fig.  76  The  best  magnetic  actuator  among  several  designed  actuators 

This  coil  actuator  ideas  was  further  tested  using  3D  model  as  shown  in  Fig.  77  . 
Considering  the  symmetry,  one  fourth  of  the  actuator  was  modeled.  The  geometry  was 
modeled  using  commercial  software.  Pro/engineering.  Fig.  78  shows  the  dimension  of 
the  3D  coil  actuator  and  the  orientation  of  the  permanent  magnets.  The  dimensions 
shown  are  all  in  mm. 


Permanent  magnets 


Fig.  77  Soid  modei  of  the  3D  coii  actuator 
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Fig.  78  The  top  view  of  the  3D  coil  actuator 

The  moving  coils  carry  current,  coils  that  has  a  conductivity  values  of  \(f  S/m.  In  order 
to  perform  the  analysis,  eight  node  brick  elements  were  used.  The  total  number  of  the 
nodes  was  9651.  According  to  the  design  criterion,  the  magnetomotive  force,  Nl,  can 
vary  from  0  to  168.  When  Nl  is  equal  to  30,  the  magnetic  flux  density  of  the  coil  actuator 
is  shown  in  Fig.  79  . 
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Fig.  79  The  magnitude  of  the  magnetic  flux  density  of  the  coil  actuator 

The  magnetic  flux  density,  B,  vary  from  1.17  to  2.33  T  nearby  the  moving  coils,  so  the 
computed  force  is  -0.114  N.  When  Nl  is  equal  to  10,  30,  60,  90,  120  or  168,  the  computed 
Lorentz  forces  are  plotted  in  Fig.  80  .  This  figures  shows  the  linear  relation  between  the 
magnetomotive  force  and  the  Lorentz  force.  The  maximum  force  is  0.683  N. 
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Fig.  80  Lorentz  force  of  the  coil  actuator  versus  Nl 

7.  Coupled  magneto-elastostatic  analysis  of  flapping  wing  actuation 

The  best  coil  actuator  is  embedded  in  a  micro  air  vehicle  with  flapping  wings  as 
shown  in  Fig.  81  .  The  coil  actuator  can  fit  inside  of  the  fuselage.  The  moving  coil  of  the 
coil  actuator  is  attached  to  the  structure  or  could  even  be  built-into  the  structure. 


Fig.  81  The  coil  actuator  with  flapping  wings 
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Four  different  structure  supports  for  the  wing  were  designed  as  shown  in  Fig.  82  .  The 
structures  are  modeled  as  shells  and  surface  models  as  shown  in  Fig.  83  used  represent 
their  geometry.  The  thickness  of  the  shell-like  structure  is  assumed  to  be  equal  to 
0.45mm.  All  the  designs  have  a  similar  top  view;  however,  they  have  different  front  and 
side  views.  The  dimensions  shown  are  all  in  mm.  As  the  length  in  chord  wise  direction 
is  26  mm,  two  coil  actuators  can  be  placed  along  the  chordwise  direction  because  the 
width  of  the  coil  actuator  is  12.2mm. 

Top  view 

H  "  160  -  - — 


D 


Fig.  82  Four  structures  with  flapping  wings.  A)  Design  1,  B)  Design  2,  C)  Design  3,  and  D) 

Design  4. 
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Fig.  83  Four  surface  structures  with  flapping  wings.  A)  Design  1,  B)  Design  2,  C)  Design  3, 

and  D)  Design  4. 

In  order  to  perform  analysis  of  thin  shell-like  structures,  IBFEM  shell  elements 
described  in  the  previous  chapter  were  used.  The  total  number  of  nodes  in  model  is  666. 
Considering  the  symmetry  of  the  geometry,  half  of  the  structure  was  modeled  for  the 
analysis.  The  structured  mesh  and  the  boundary  conditions  are  shown  in  Fig.  84  .  The 
magnetic  force  acts  downward  so  that  the  wing  undergoes  an  upstroke.  When  one  coil 
actuator  is  used,  the  applied  force  varies  from  0  to  0.342  N  at  the  edge.  Using  two  coil 
actuators,  the  magnetic  force  can  be  double. 
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Natural  boundary  condition 


Essential  boundary  condition 


Fig.  84  Structured  mesh  and  boundary  conditions  of  the  first  design 

By  reversing  current  direction,  force  can  be  reversed  so  that  the  wing  can  undergo  a 
down-stroke.  Using  the  coil  actuator,  the  same  amount  of  Lorentz  force  can  be 
produced  in  the  opposite  direction  by  changing  the  direction  of  the  current.  When  NI  is 
equal  to  30,  the  wing  up-strokes  and  down-strokes  of  the  four  designs  are  shown  in  Fig. 
85  .  For  the  up-stroke,  the  maximum  displacement  on  the  tip  is  2.633  xl0"^mw  in  the 
first  design,  4.238 x  10“^ ww  in  the  second  design,  3.765 x  10“^ ww  in  the  third  design,  and 
1.968xlO“'mOT  in  the  fourth  design.  Magnitudes  of  the  tip  deflections  during  the  wing 
up-stroke  are  the  same  as  ones  during  the  wing  down-stroke  because  the  magnetic  force 
is  proportional  to  the  displacement.  Thus,  the  tip  deflection  can  be  double  including 
both  strokes.  Among  the  four  designs  studied  here,  the  last  design  produces  the  largest 
tip  deflection. 
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Fig.  85  Displacement  in  the  z-direction  during  the  wing  stroke.  A)  Design  1,  B)  Design  2,  C) 

Design  3,  and  D)  Design  4. 
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Tip  displacement  [mm] 


varies  from  1  to  168,  the  tip  deflection  as  a  function  of  the 
own  in  Fig.  85  .  The  first  three  designs  are  too  stiff  for  our  app 


VII.  CONCLUSIONS  AND  FUTURE  RESEARCH  DIRECTION 


1.  Summary 

The  goal  of  this  project  was  to  explore  means  of  designing  multi-functional 
composite  structures  that  have  sensing  and  actuation  capabilities.  Both  sensing  and 
actuation  means  were  studied.  Computational  tools  were  created  that  help  in  designing 
and  analyzing  magnetically  actuated  composite  structures.  A  brief  summary  of  each  of 
these  activities  is  given  below. 

Composite  sensors: 

A  method  for  detecting  load  application  and  damage  in  composite  structures  by 
determining  change  in  resistivity/conductivity  was  studied.  An  algorithm  for  solving 
inverse  problem  to  determine  average  resistivity  values  in  composite  structures  was 
demonstrated.  The  method  can  use  data  from  an  arbitrarily  large  number  of  electrodes 
to  compute  average  values  of  resistivity  or  conductivity  for  the  structure.  Finite  element 
models  for  the  structure  are  used  to  solve  the  forward  problem,  making  this  method 
very  general  and  applicable  to  arbitrary  shaped  structures.  Ideally,  the  electrodes 
should  be  embedded  in  the  structures  during  the  manufacturing  process  itself  so  that  it 
can  be  used  for  quality  control,  detection  of  defects  as  well  as  subsequent  health 
monitoring.  One  of  the  advantages  of  measuring  resistivity  is  that  damage  can  be 
detected  even  in  structures  that  were  not  tested  during  manufacturing.  Damage  can  be 
detected  for  structures  that  are  in  use  by  attaching  electrodes  on  the  surface, 
determining  the  average  resistivity  and  comparing  it  to  values  associated  with 
undamaged  material.  The  main  source  of  error  in  this  approach  arises  from  inaccuracy 
in  the  geometric  models  of  the  structure  and  the  electrodes.  Random  noise  added  to  the 
voltage  data  used  in  numerical  validation  indicates  that  any  error  in  the  voltage  data 
can  get  amplified  due  to  difficulties  in  numerical  convergence.  In  principle,  a  similar 
inverse  approach  can  also  be  used  to  determine  applied  loads  on  the  structure. 
Experimental  studies  indicate  that  this  may  be  difficult  because  the  changes  in  the 
electrode  voltage  due  to  strains  are  very  small.  Even  with  amplification,  the  data  is  hard 
to  use  because  of  significant  non-linearity  in  the  observed  behavior.  However,  the 
approach  we  developed  is  promising  for  detecting  damages  or  defects  because  they 
cause  significant  changes  in  resistivity  and  is  therefore  easier  to  detect.  Further  study  is 
needed  to  explore  ways  of  determining  applied  loads  and  strains. 
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Magnetically  actuated  composite  structures 

The  idea  of  using  magnetic  forces  to  actuate  structural  mechanisms  was  studied. 
The  main  application  of  interest  is  micro  air  vehicle  wings  that  are  shell  like  structures. 
Topology  optimization  method  was  studied  as  a  potential  method  for  designing 
structures  that  have  specified  modes  of  deformation.  The  structure  is  then  to  actuated 
using  magnetic  actuation  means  built  into  or  around  the  structures.  Several  actuators 
were  studied  including  solenoid  actuators  and  coil  actuators.  After  systematic 
comparison  of  several  designs,  it  was  concluded  that  a  coil  actuator  built  into  composite 
structures  is  an  ideal  means  for  actuation  of  composite  structures.  A  conceptual  design 
of  a  flapping  wing  vehicle  was  developed  that  is  designed  to  actuate  by  built  in 
actuation  capability  of  the  body,  wing  and  support  structures.  No  external  mechanisms, 
motors  or  linkages  are  needed. 

Computational  tools  for  designing  magnetically  actuated  structures 

Computational  tools  were  developed  to  design  and  analyze  structures  actuated  by 
magnetic  forces.  Magnetostatic  analysis  capability  was  implemented  into  a  pre-existing 
software  (named  IBFEM)  developed  at  the  University  of  Florida  that  can  perform  finite 
element  analysis  without  the  need  for  generating  mesh.  Solid  and  surface  geometry 
modeled  on  commercial  CAD  software  can  be  imported  into  this  software  and  analysis 
can  be  performed  without  approximating  the  geometry  using  a  conforming  mesh.  The 
structured  mesh  approach  has  been  demonstrated  to  work  for  magnetostatic  analysis 
and  validated  using  several  examples  with  known  solutions.  The  approach  has  been 
demonstrated  for  both  2D  and  3D  magnetostatic  models.  Structured  mesh  is  easy  to 
generate  and  the  elements  are  regular  and  not  distorted  as  in  traditional  finite  element 
mesh.  Magnetic  forces  were  computed  by  integrating  the  magnetic  force  density.  These 
forces  and  then  used  in  a  subsequent  structural  analysis  to  determine  the  deflection  of 
the  structure.  Shell  elements  based  on  uniform  B-spline  shape  function  were 
implemented  into  IBFEM.  One  of  the  key  advantages  of  using  these  elements  is  that  a 
structured  mesh,  which  is  easy  generate  automatically,  can  be  used  for  the  analysis. 
Both  quadratic  and  cubic  B-spline  shape  functions  were  tested  and  it  was  found  that 
cubic  elements  provided  very  good  results  with  fewer  elements  than  quadratic 
elements.  Computational  cost  is  higher  for  these  elements  compared  to  traditional  shell 
elements  but  often  fewer  elements  are  needed  to  get  accurate  results  with  cubic 
elements.  The  time  taken  to  create  the  model  is  significantly  lower  because  structured 
mesh  generation  is  easily  automated. 


2.  Future  work 

Design  methodology  developed  here  based  on  topology  optimization  can  be 
extended  to  design  structures  that  morph  or  oscillate  in  specified  modes.  This  would 
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allow  us  to  design  wings  that  are  designed  to  oscillate  in  a  mode  shape  that  produces 
the  most  lift  and  thrust.  High  speed  photography  has  been  used  to  study  motion  of 
insect  and  bird  wings  by  biologists.  It  is  challenging  to  design  wings  that  can  oscillate  in 
a  fashion  that  mimics  the  flapping  of  insect  or  bird  wings.  Preliminary  results  obtained 
using  topology  optimization  in  this  research  suggests  that  this  is  a  promising  area  for 
further  study  but  it  is  a  challenging  problem  that  may  take  2-3  years  of  work  especially 
if  we  also  want  to  create  physical  models  for  validation.  In  order  to  fully  understand  the 
interaction  of  flapping  wings  with  air,  fluid  structure  interaction  model  are  needed. 
Incorporating  CFD  into  IBFEM  software  would  make  simulation  and  analysis  of  such 
systems  much  easier  than  currently  available  tools. 

The  magnetostatic  analysis  capability  that  was  developed  in  this  research  needs  to 
be  extended  in  the  future  to  nonlinear  problems  as  well  as  dynamics  problems.  Since  3D 
analysis  is  computationally  expensive,  it  is  not  always  desirable  to  use  a  uniform  mesh 
everywhere.  Mesh  refinement  techniques  are  needed  that  locally  refines  the  mesh  while 
still  using  a  structured  mesh  that  has  regular  shaped  elements.  Development  of  such 
local  refinement  techniques  is  part  of  the  future  work  needed  to  make  3D  magnetostatic 
analysis  less  expensive.  For  application  to  structural  mechanisms  undergoing  large 
deformation,  the  shell  elements  must  be  extended  to  allow  large  elastic  deformation  and 
elasto-plastic  deformations.  The  performance  of  the  IBFEM  shell  elements  can  also  be 
further  improved  by  incorporating  adaptive  mesh  refinement.  The  shell  elements  using 
IBFEM  can  also  be  extended  to  mixed  solid-shell  type  analysis  where  thin  shells  are 
attached  to  solid  structures. 

This  research  focused  on  the  steady  static  analysis  so  that  the  coupled  problem  is 
modeled  as  a  weakly  coupled  magneto-elastostatic  analysis.  This  assumption  allows  us 
to  perform  the  magnetostatic  and  elastostatic  analysis  sequentially.  If  the  magnetic  field 
changes  due  to  the  elastic  deformation,  the  problem  becomes  nonlinear  and  a  strongly 
coupled  nonlinear  analysis  is  needed.  The  magnetic  field  may  change  because  the 
structural  deformation  may  cause  permanent  magnets  or  circuits  attached  to  the 
structure  to  also  move  relative  to  each  other.  If  the  structure  undergoes  large 
deformation  then  geometric  nonlinearities  must  be  included.  In  addition,  if  the  coupled 
problem  is  a  dynamic  problem,  then  the  electrical  circuits,  the  magnetic  circuits  and  the 
structural  analysis  must  be  solved  as  a  coupled  problem.  Another  source  of  nonlinearity 
is  due  to  contact  between  the  armature  and  the  stator. 
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